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ABSTRACT 

This  report  investigates  the  stability  properties  of 
a  screw  pinch  equilibrium  with  flow.   The  question  of 
local  stability  is  fully  answered  by  the  derivation  of  a 
necessary  and  sufficient  criterion.   Unlike  the  static 
case,  local  modes  can  be  suppressed  by  global  properties 
of  the  plasma,  as  well  as  by  a  large  sheared  flow.   When 
local  instability  does  occur,  the  accumulation  point  of 
the  unstable  modes  is  found  to  be  the  point  of  overlap 
of  the  Alfven  and  cusp  continua. 

The  last  chapter  of  the  work  develops  a  sufficient 
(but  not  necessary)  condition  for  stability.   Some  bounds 
on  the  growth  rates  of  instabilities  are  also  obtained. 
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Introduction 

The  term  "stability"  in  the  context  of 
magnetohydrodynamic  systems  is  usually  used  to  describe 
a  situation  where  small  perturbation  of  a  given  equilibrium 
state,  does  not  result  in  large  change  of  the  system.   A 
simplified  version  of  this  concept  is  linear  stability. 
It  is  assumed  that  the  disturbances  are  so  small  that  a 
linearized  version  of  the  equations  of  motion  is  suf- 
ficient to  describe  the  system,  and  for  stability  it  is 
required  that  solutions  of  the  linearized  equations  do 
not  grow  unboundedly  in  time.   Mathematical  analysis  is 
simplified  if  we  ignore  modes  that  do  not  grow  exponen- 
tially in  time.   We  adopt  this  common  practice  here. 

The  model  we  use  in  this  thesis  is  ideal  MHD,  in 
which  plasma  is  described  as  a  non  dissipative  fluid, 
acting  under  the   influence  of  hydrodynamic  and  electro- 
magnetic forces.   We  use  one  of  the  simplest  configurations 
for  containment  of  plasma  -  a  diffuse  screw  pinch,  in 
which  the  plasma  occupies  the  volume  of  a  cylinder  and  is 
confined  by  the  current  flowing  within  it.   Our  model 
allows  for  an  equilibrium  state  with  steady  flow,  and  the 
investigation  of  this  state  is  what  distinguishes  this 
work  from  others  previously  done.   A  quick  review  of  these 
works  would  provide  us  with  a  deeper  insight  into  the 
difficulties  involved. 


The  question  of  stability  of  MHD  equilibria  without 
flow  was  answered  in  principle  by  Bernstein  et  al  [1] . 
These  authors  used  an  energy  principle  which  provided  a 
necessary  and  sufficient  condition  for  exponential 
stability.   The  necessity  part  was  proven  by  Laval, 
Mercier  and  Pellat  [20].   The  energy  principle  was  worked 
out  in  a  straightforward  way  by  Newcomb  [11]  to  provide 
a  marginal  stability  criterion  for  the  static  screw  pinch 
equilibrium.   Goedbloed  [21]  arrived  at  similar  results 
using  a  different  argument  based  on  oscillation  theorems. 
The  success  in  solving  the  static  problem  was  made  pos- 
sible by  the  availability  of  a  self  adjoint  formulation 
of  the  problem.   The  general  solution  to  the  equation 
u    +  HU  =  0,  where  u  belongs  to  some  Hilbert  space  and 

H  is  a  self  adjoint  operator  in  that  space,  is  given  by 

*.u   i  •         i- •  i/H  t   ,,  -i/H  t    , 

the  linear  combination  u  =  ae     u.+be      u_  where 

the  exponential  operators  are  defined  as  the  appropriate 

integral  taken  over  the  spectrum  of  H  and  u1  ,  u-  are 

vectors  in  the  space.   Clearly  H  must  be  positive  for 

exponential  stability.   This  demonstration  will  convince 

the  reader  of  the  usefulness  of  the  self  adjoint 

formulation  (although  H  in  reference  1  was  only  shown  to 

be  formally  self  adjoint) . 

When  flow  is  introduced  into  the  system,  the  self 

adjointness  is  lost  and  there  is  a  jump  in  the  degree  of 

difficulty  in  treating  the  case.   Only  a  few  steady 


equilibria  were  investigated  prior  to  this  time. 
Chandrasekhar  [22]  and  Newcomb  [23]  found  stability 
criteria  for  purely  rotational  flows  with  a  uni- 
directional magnetic  field,  and  with  respect  to  azimu- 
thally  symmetric  perturbations  only.   Spies  [24] 
recently  worked  out  the  stability  of  the  6-pinch 
(axial  magnetic  field)  with  solid  rotational  flow 
and  with  respect  to  axially  independent  perturbations. 

In  this  work  we  consider  equilibria  with  magnetic 
fields  having  toroidal  as  well  as  poloidal  components, 
and  flows  having  axial  and  azimuthal  components.   The 
price  we  pay  for  this  generality  is  in  being  unable  to 
conclude  categorically  whether  a  system  is  stable  or 
not,  but  we  have  to  focus  our  attention  on  parts  of 
the  spectrum. 

In  chapter  II  we  investigate  the  continuous 
spectrum  corresponding  to  the  problem,  and  find  that 
it  consists  of  real  members  only  and  thus  does  not 
cause  exponential  instability.   Another  part  of  the 
spectrum  that  is  relatively  easy  to  handle  is  the 
accumulating  sequences  of  eigenmodes,  known  as  local 
(Suydam)  modes.   The  investigation  of  such  unstable 
modes  requires  the  development  of  an  asymptotic  ex- 
pansion technique  for  a  non-standard  case  of  a  singular 
perturbation.   This  is  carried  out  in  chapter  III  and 


appendix  A.   It  is  then  found  in  chapter  IV  that,  like 
the  static  case,  unstable  eigenmodes  can  only  accumu- 
late at  the  point  of  overlap  of  the  various  continuous 
spectra.   We  derive  a  necessary  and  sufficient  condition 
for  local  stability,  which  unlike  Suydam's  condition 
for  the  static  case,  is  not  entirely  local  and  depends 
on  global  properties  of  the  plasma  which  may  suppress 
the  local  modes.   Another  stabilizing  effect  is  generated 
by  the  shear  of  the  equilibrium  flow.   A  large  enough 
flow  shear  is  sufficient  for  local  stability. 

In  chapter  V  a  Lagrangian  description  of  the  motion 
is  used  to  obtain  some  bounds  on  the  growth  rates  of 
unstable  modes  and  we  develop  a  sufficient  condition 
for  stability,  which  reduces  to  Newcomb ' s  [11]  necessary 
and  sufficient  condition  in  the  static  case. 


CHAPTER  I:   Basic  Equations 

1.1   The  Eulerian  Equations  of  Motion 

The  equations  which  describe  an  ideal  MHD  (i.e. 
perfectly  conducting)  plasma  are: 


ft  +  P  div  u  =  0 


(MF) 


du 
p  —  =  -  grad  p+j_xB 


dt 


§f-° 


3B 

^—  +  curl  E  =  0 


div  B  =  0 


curl  B  =  2. 


E+uxB  =  0 


p  =  f(p,S) 


where  p,  p,  S,  u,  E,  B,  j_  are  the  plasma  mass  density, 
pressure,  entropy,  velocity,  the  electric  and  magnetic 
fields  and  the  current  density.   The  symbol  -T-r-  stands 
for  the  Lagrangian  time  derivative  —  +  u  •  '. '  . 

d  t     —    — 


The  first  three  equations  express  the  conservation  of 
mass,  momentum  and  entropy.   The  next  three  equations  are 
Maxwell's  equations  for  the  electromagnetic  field  with  the 
displacement  current  omitted.   The  last  two  equations  are 
Ohm's  law  for  a  perfectly  conducting  fluid   and  the  given 
equation  of  state.   The  conditions  for  the  applicability 
of  these  equations  are  discussed  in  reference  1. 

We  use  the  last  three  equations  to  eliminate  E,j_  and 
S  from  the  system.   In  particular,  we  multiply  the  first 
equation  by  |£  and  the  third  one  by  ^    •   Adding  these 
two  equations  yields  an  equation  for  p  which  replaces  the 
adiabatic  law.   System  (MF)  takes  now  the  simplified  form 
of  the  Lundquist  equations  [2]  : 


U  +  p  div  u  =  0 
dt         — 


(L) 


-£■   +  pa   div  u  =  0 
dt  — 


du  ,  2 

P  3t  =  "  9rad(p  +  jB  )  +  (B-V)B 

9B 

—  +  curl(Bxu)  =  0 

at      — 

div  B  =  0 


where  a   =  |^  is  the  square  of  the  speed  of  sound.   It  is 
well  known  that  the  last  equation  can  be  considered  an  initial 


condition,  since  it  will  be  satisfied  for  all  time  if 
satisfied  initially.   We  point  out  that  system (L)  is 
hyperbolic  [3] ,  a  useful  property  for  our  work. 

In  this  treatise  we  are  interested  in  a  special 
plasma  configuration:  the  diffuse  screw-pinch  with  flow. 
The  following  assumptions  characterize  the  steady  state: 
i.      The  plasma  occupies  the  volume  of  a  cylinder 

of  length  L  (along  the  z  axis,  say)  surrounded 
by  a  perfectly  conducting  wall, 
ii.     All  quantities  (the  pressure,  magnetic  field, 
etc.)  depend  on  the  distance  from  the  axis 
only  (r  -  in  cylindrical  coordinates) . 
iii.    The  equilibrium  magnetic  field,  and  consis- 
tently the  flow  velocity,  have  vanishing 
radial  components. 
Such  a  flow  has  to  satisfy  only  one  non  trivial 
equation  of  the  system  (L) 

B2     u2 

(i)  (P  +  |b2)'+  4=  p^| 

where  the  prime  denotes  differentiation  with  respect  to  r 
and  subscripts  denote  vector  components  in  cylindrical 
coordinates.   Since  the  axial  velocity  does  not  enter  this 
equation,  one  equilibrium  state  gives  rise  to  many  others 
by  simply  changing  u  . 


The  boundary  conditions  that  have  to  be  imposed 
consistently  with  the  cylindrical  geometry,  are 

i.      Periodicity  in  z  with  a  period  L. 

(BC)  ii.     B   =  u   =  0  at  the  outer  wall, 
r     z 

iii.    Regularity  at  the  origin. 
These  conditions  apply  of  course  to  steady,  as  well  as 
to  time-dependent  flows.   The  second  condition  is 
automatically  satisfied  by  the  steady  state. 


1.2   The  Perturbed  Equations 

We  turn  now  to  the  stability  analysis  of  the 
system  (L) .   Rather  than  working  with  the  non-linear 
system,  we  restrict  our  attention  to  the  question  of 
linear  stability  of  this  system.   Accordingly  we  consider 
a  small  perturbation  of  the  equilibrium  state,  and 
linearize  the  system(L)  about  it.   The  coefficients  of 
the  linearized  equations  then  depend  on  r  only,  and  it 
is  possible  to  Fourier  transform  the  system  with  respect 
to  the  axial  and  angular  variables,  and  to  Laplace 
transform  it  in  time.   This  is  equivalent  to  assuming  a 
dependence  of  every  perturbed  quantity  E, ,    of  the  form 

|(r,6,z,t)  =  5(r)  exp  i(m8+kz+ut). 

Application  of  the  Laplace  transform  is  justified  by  the 
hyperbolocity  of  the  system  (L) ,  from  which  we  conclude 
that  the  linearized  system  is  of  a  similar  type  and  thus 
possesses  only  solutions  that  grow  at  most  exponentially 
fast  in  time. 

The  distinction  between  the  Laplace  and  Fourier 
transforms  is  made  in  order  to  underline  the  existence  of 
non-homogeneous  terms  originating  from  initial  conditions. 
These  terms  play  an  important  role  in  identifying  the 
continuous  spectrum.   However,  since  our  main  interest  is 
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the  point  spectrum,  we  actually  omit  such  terms. 

On  linearization  about  the  screw  pinch  with  flow  as 
the  steady  state,  the  system  (L)  can  be  reduced  to  a  2x2 
system  of  first  order  ordinary  differential  equations  in 
only  two  of  the  perturbed  quantities.   All  other  per- 
turbed quantities  can  be  expressed  in  terms  of  these  two. 
This  substantial  simplification  is  made  possible  by  the 
choice  of  the  screw  pinch  as  our  equilibrium  state,  which 
in  turn  makes  the  linearized  system  (L)  algebraic  (and 
not  differential)  in  six  of  the  eight  dependent  variables. 
The  reduced  system  can  be  written  as 


(A')  aAC("   )'  +  A°(^   )  =  0 


where  u=u   is  the  radial  component  of  the  perturbed  flow, 
and  p*  =  p+B  •  B  is  the  perturbed  total  pressure.   All 
other  perturbed  quantities  can  be  expressed  as  linear 
combinations  of  u,  u1,  p*.   The  explicit  dependence  is 
described  in  Appendix  B. 

The  coefficients  of  system  (A1)  are  functions  of  the 
unperturbed  flow.   To  avoid  confusion,  an  equilibrium 
quantity  will  carry  a  subscript  o  when  necessary.   We  use 
the  following  definitions: 

o   =  u+k'u 

2    2 

A  =  H  -pa* 

2    2    2 
C  =  B*(  H  -pa  ) 
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and  we  use  the  notation 


k  =  (0,— ,k)  in  cylindrical  coordinates 


K2  =  k2  +  Hi! 

r 


and         u  =  (0,v,w)  is  the  equilibrium  flow, 


H  =  k-B 


2      2   2 

B*  =  pa^+B^ 


2 
j 

b;     a"+b' 


r  -  Pa       a 
~2 T~2 


2         B2 
and       b   =  —   is  the  Alfven  speed , 


o         1  n2 


o    2   o 


The   matrix   elements   of    A      are 


Xll    =    -C[^T   a(BeH-pov)+Ao'-    p-Ao]-   pp2o3v2A   +    pp2o4Be  (Hv-BQo) 

X12    =    °2<P2°4    +    CI<2) 

A°      =      C[Ar(|   p°')'    +    A2    -    -j(B    H    -    pov)2]     -    [  ±-   pv2A 

r 

-    -   poBfl(Hv-B    J))2 
r  0  0 

A°       =    -A°       -    AC 0 '     +    -ACo 
22  11  r 
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System  (A)  is  supplemented  with  two  boundary  conditions 

(Al)     u  =  0  on  the  wall  r=  R 

(A2)     ru  is  regular  near  the  origin  r  =  0. 

Every  system  of  the  form  (A)  can  be  reduced  to  a 
single  second  order  equation  in  any  of  the  unknowns,  e.g 


(A)  (pu') '  +  qu  =  0 


rAC 
where  p  =  

A° 
A12 


o  A . 


q  .  £[  det  A     +   (  -11).] 
aAC.A12         A12 


Although  equation  (A)  has  a  self  adjoint  form,  it  does 

not  have  this  property.   The  coefficients  p,q  contain 

03  which  may  be  complex,  thus  destroying  the  self  adjointness 


with  respect  to  the  relevant  scalar  product  (u,v)  =    uvdr 

It  should  be  pointed  out  that  every  second  order  linear  ODE 
can  be  brought  into  form  (A)  by  multiplying  with  an  appropri- 
ate factor,  and  that  the  self  adjoint  form  has  no  meaning 
of  its  own.   We  choose  it  for  reasons  of  convenience  in 
handling  the  associated  Green's  function,  and  also  in  order 
to  show  the  similarity  to  the  equation  derived  by  Hain 
and  Lust  [4]  for  the  static  screw  pinch.   The  two  equations 
are  in  fact  identical  when  we  let  the  equilibrium  flow 
vanish. 
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CHAPTER  II:   The  Spectrum 

II. 1   Definitions 

The  term  "spectrum"  is  properly  used  in  the  context 
of  operators  in  a  topological  space  -  most  commonly 
Hilbert  space.   The  present  formulation  of  the  stability 
problem  deals  with  the  growth  in  time  of  solutions  to  the 
linearized  Lundquist  equation.   We  did  not  introduce  a 
space  to  which  the  eight  initial  perturbations  p,p,u,B 
had  to  be  restricted,  and  in  fact  there  does  not  seem  to 
be  a  "natural"  space  associated  with  this  problem,  with 
a  norm  that  would  for  example  provide  a  spectral  expansion 
or  have  other  advantages  of  that  kind.  The  energy  norm 

t 

p  |u|  dV  which  is  used  in  the  static  case,  would  limit 
consideration  to  perturbations  of  finite  kinetic  energy 
but  would  offer  no  mathematically  exploitable  structure 
to  the  problem  at  hand.   The  only  clue  we  have  is  that 
the  "right"  norm,  if  there  is  any,  must  measure  the  growth 
of  three  perturbed  quantities,  since  they  all  depend  on 
pA,u,u'.   It  will  be  shown  in  chapter  V  that  the  Lagrangian 
formulation  of  the  problem  does  offer  an  elegant  way  of 
turning  the  stability  question  into  a  non-linear  eigenvalue 
problem  in  a  Hilbert  space.   In  the  meantime,  however,  we 
restrict  the  discussion  to  a  purely  differential-equation 
context,  although  we  continue  to  use  the  term  "spectrum" 
for  simplicity  and  after  proper  definition. 
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Our  guide  for  the  use  of  this  term  is  the  Laplace 
transform  which  was  applied  in  the  normal  mode  analysis, 
Every  perturbed  quantity  £(r,to)  which  is  a  solution  of 
the  non-homogeneous  Laplace  transformed  equation,  cor- 
responds to  a  time  dependent  perturbation  £ ,  obtainable 
by  the  inverse  transform 


«<r'fc>  =  h 


i  cot  „  ,    ,  , 
e    £  (r  ,  to)dto   , 


r 


where  T  is  a  horizontal  line  in  the  complex  10  plane 
lying  below  all  the  singularities  of  £(r,to).   The  hyper- 
bolicity  of  the  problem  guarantees  that  such  a  line  T 
exists,  i.e.  that  £(r,to)  is  analytic  in  to  for  Im  to  sufficient- 
ly negative.   We  can  now  continue  £(r,to)  analytically  and 
move  the  line  of  integration  upward.   As  we  pass  a 
singularity  of  E, ,    being  either  an  isolated  singularity 
or  a  singularity  along  a  curve,  we  have  to  add  to  the 
integral  along  T,  an  integral  around  the  singularity  we 
passed.   This  process  yields  a  "spectral  resolution"  of 
the  time  dependent  solution  in  the  eigenmodes  of  the 
system  of  equations,  and  accordingly  we  call  the  isolated 
singularities  of  £  -  "eigenvalues"  of  the  problem,  and  the 
other  singularities  -  the  "continuous  spectrum".   The  rest 
of  this  chapter  is  devoted  to  the  identification  and  in- 
vestigation of  the  different  parts  of  the  spectrum  as- 
sociated with  the  screw  pinch  problem. 
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II. 2   The  Continuous  Spectrum 

In  all  subsequent  discussions  we  consider  the  wave 
numbers  m  and  k  to  be  fixed.   System  (A')  indicates  the 
existence  of  a  continuous  spectrum  which  contains  all 
points  oj  for  which  the  system  is  singular  at  some  point 
r  in  the  range  of  the  plasma.   Although  the  spectrum  is 
associatied  with  solutions  of  the  non-homogeneous  prob- 
lem, it  can  be  detected  from  the  non-analyticity  of 
solutions  of  the  homogeneous  problem  simply  because  the 
Green's  function  of  the  problem  is  constructed  by 
solutions  of  the  homogeneous  equation.   There  are  three 

parts  to  the  continuous  spectrum: 

2    2 
i.       The  Alfven  continuum:       A  =  H  -pa   =  0 

2     2 
ii.      The  cusp  continuum:        O  8H   -pa   =  0 

iii.     The  Eulerian  continuum:  a  =  0. 

It  should  be  pointed  out  that  although  points  where  X    _  =  0 

appear  to  be  singularities  of  equation  (A) ,  they  are  in 

fact  regular  points,  as  seen  from  the  equivalent  system  (A1) 

This  apparent  singularity  results  from  the  elimination  of 

p^  as  an  unknown.   If  we  eliminate  u  instead,  we  will  get 

Xj-x    as  an  apparent  singularity. 

The  first  two  continua  are  known  in  the  static 

equilibrium  case [5],  and  have  there  identical  expressions, 

with  o  replaced  by  w.   The  last  continuum  reduces  to  the 

point  uj=0  when  the  flow  vanishes,  and  is  also  a  known 

singularity  in  the  static  case.   The  source  of  the  names 
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assigned  to  the  first  two  continua  is  their  similarity 
to  characteristic  speeds  of  the  Lundquist  equation  [3]. 
The  Alfven  continuum  is  defined  by  an  expression 
identical  to  the  Alfven  wave  propagation  speed  in  the 
direction  k.   The  cusp  continuum,  if_  treated  as  a 
propagating  signal,  would  yield  propagation  along  a 
magnetic  field  line  with  a  speed  equal  to  that  of  the 
trailing  cusp  of  the  slow  magnetosonic  wave.   It  will 
be  shown  in  chapter  V  that  the  Eulerian  continuum  is 
characteristic  of  the  Eulerian  description  of  the  motion 
and  is  missing  when  Lagrangian  formulation  is  used. 

Because  of  the  apparent  similarity  between  continua 
and  wave  propagation  speeds,  it  should  be  stressed  that 
these  two  phenomena  are  not  directly  related,  as  pointed 
out  by  Grad  [5] .   Firstly,  the  fast  and  slow  waves  do  not 
have  counterparts  among  the  continua.   Secondly,  the  cusp 
continuum  is  not  a  signal  speed.   The  resemblance  can  be 
explained  by  the  formal  mathematical  procedures  that  lead 
to  the  determination  of  characteristic  speeds  and  continua, 
In  deriving  the  characteristics  of  a  first  order  partial 
differential  equation,  we  make  the  change  •* —  -*■  in  6u  for 
every  differential  term,  with  n   a  constant.   This  looks 
much  like  normal  mode  analysis  which  proceeds  similarly, 
except  that  r-derivatives  are  kept,  and  non-differentiated 
terms  are  not  dropped.   The  continuum  is  determined  by  the 
vanishing  coefficient  of  the  r-derivative .   Thus,  for  w  in 
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the  continuum  much  of  the  difference  between  the  two 
processes  is  eliminated  because  there  are  no  r-differen- 
tiated  terms,  and  hence  the  similarity. 

Before  concluding  this  part,  we  observe  that  the 
continuous  spectrum  consists  only  of  real  numbers  w  and 
therefore  does  not  cause  exponential  growth  or  decay  of 
the  perturbation.   The  question  of  exponential  stability 
is  thus  decided  by  the  point  spectrum,  and  we  proceed 
with  its  investigation. 
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II. 3   Accumulating  Eigenmodes 

We  first  make  the  useful  transformation 
(2)   u  =  ia£ 

which  uses  £  ,  the  Lagrangian  radial  displacement,  as  the 
new  unknown.   The  transformation  fails  when  o  vanishes, 
i.e.  when  co  is  in  the  Eulerian  spectrum,  but  the  point 
spectrum,  which  is  the  subject  of  this  investigation,  is 
not  affected.   System  (A1)  now  takes  the  form 


(B«)  AC(^  )  '  +  A(|  )  =  0 


where  the  matrix  elements  are 


All  =  -C[^(B0H-pav)-  h]    +   |p2a3BQ  (Hv-BQo)  -  ±-  p2a2v2A 

2  4    2 
X12  =  P  °  +CK 


X2]_  =  C[Ar(i  pj')  '+A2 i(BeH-pov)2]-[^Pv2A-  ^paBg  (H v-BqO)  ] 

X    =  -A..  +  -AC 
22      11    r 


or  equivalently 
(B)  (pC)  '  +  q£  =  0 


«  -  rAC 

P  "  x — 

A12 

=  r[  det  A  +  (_li)  •] 

ACX12     X12 


2 
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The  boundary  conditions  become 

(Bl)   5  =  0  at  the  outer  wall  r  =  R. 

(B2)  rE,    is  regular  near  the  origin. 

The  transformation  to  the  Lagrangian  variable 
conceals  the  spectrum  o  =  0  and  makes  equation  (B) 
easier  to  handle.   It  does  not  however  eliminate  this 
continuum  from  the  system.   When  solving  the  initial 
value  problem  we  still  have  —  terms  in  the  non  homo- 
geneous side  of  (B) .   We  also  have  such  terms  when 
expressing  the  other  perturbed  quantities  (B,  etc.) 
in  terms  of  pt,  E,   and  £ '  . 

The  point  spectrum  consists  of  all  numbers  oj  not 
in  the  continuum,  for  which  there  exists  a  non  trivial 
solution  to  equation  (B) ,  which  satisfies  the  boundary 
conditions  (B1-B2).   Clearly,  any  non  real  contribution 
to  the  coefficients  p,q  of  (B)  is  due  to  u>.   By  taking 
the  complex  conjugate  of  (B)  we  find  that  if  u  is  an 
eigenvalue,  so  io  its  complex  conjugate  u>  (with  an 
eigenf unction  \  )  . 
Corollary :   The  flow  is  exponentially  stable  if  and  only 

if  there  are  no  non-real  eigenvalue*. 

Rather  than  investigating  the  question  of  the  existence 
or  non-existence  of  complex  eigenvalues,  we  tend  to  the 
simpler  problem  of  whether  or  not  there  exists  an  accumulat- 
ing sequence  of  unstable  eigenva  .istence  of 
such  a  sequence  will  be  referred  to  as  "lo< 
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for  reasons  that  will  be  made  clear  in  chapter  III.  Our 
interest  in  local  stability  stems  from  the  fact  that  in 
the  steady  screw  pinch  case  local  stability  indicates 
global  stability  in  current  experimental  ranges  of  the 
plasma  parameters  [6] . 

One  way  of  finding  the  eigenvalues,  is  to  start  with 
a  solution  £  of  (B)  which  satisfies  (B2)  at  r  =  0 ,  and  in 
addition  satisfies  another  arbitrary  condition  at  the 
origin  independent    of  w.   The  eigenvalues  are  then  the 
roots  of  the  equation  £(R,co)  =  0,  which  is    condition 
(Bl)  .   £(R,co)  is  an  analytic  function  of  to  for  go  not  in 
the  Alfven  or  cusp  continuum.   The  principle  of  isolated 
zeroes  of  analytic  functions  implies  that  the  zeroes  can 
only  accumulate  at  a  singularity  of  the  function,  namely 
the  Alfven  and  cusp  continua,  or  possibly  the  point  at 
infinity. 

Claim:   The  point  at  infinity  is  not  an  accumulation 
point  of  complex  eigenvalues . 

Proof:   If  a)  is  an  eigenvalue  with  an  eigenf unction  £(r,io), 
we  multiply  equation  (B)  by  f  and  integrate  by  parts. 
This  yields 


(3) 


R 

«  I  2  _  „l  r  I  2' 


(pU '  I     -  q|sr>dr  =  o. 


There  is  no  contribution  from  the  boundary  because  £ 
satisfies  (B1-B2) .   The  coefficients  p,q  of  (B)  are 
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rational  functions  of  a)  and  possess  the  following 
expansion  near  a)  =  °°: 


(4) 


P  =  Po+P2  4    +  0(-T) 

0)  0) 


q  =  q  w   +  0(a)) 
o 


where  p  ,p0,q   are  positive  functions  for  all  r  4=  0 
o  2.      o     r 

and  are  given  by 


P0  -  rB* 
p2  =  -  B^(K2B^  -  (3+l)H2) 


qo  =  rp 


2   2 
(We  add  the  assumption  that  p,  a  ,B*  do  not  vanish 

anywhere  in  the  plasma).   The  imaginary  part  of  (3) 

has  the  form 


(5) 


(Im  p)  \C   |2  "  dm  q)  |£|2  =  0 


For  |  Co  |  large  enough  but  not  real,  Im  q  and  -Im  p  have 
the  same  sign  except  when  a)  is  near  the  imaginary  axis. 
Therefore  (5)  can  only  be  satisfied  for  such  oj.   To  ex- 
clude this  possibility  we  take  the  real  part  of   (3) . 

(6) 


(Re  p)  |£'  |2  -  (Re  q)  \E,\2    =  0 


Re  q  is  negative  for  uj  large  enough  and  near  the  imaginary 
axis,  while  Re  p  is  positive  and  (6)  is  impossible. 
The  last  proof  excludes  accumulation  near 

infinity  of  complex  eigenvalues  only.   It  does  not 


22 


exclude  accumulation  of  real  eigenvalues.  Indeed,  the 
asymptotic  behavior  of  solutions  of  (B)  for  large  w  is 
[7] 


£  =  cos  [co 


r 

q  1  . 

(-^;  2  dr]  [f  (r)+0(±)] 
PQ  w 


which  verifies  the  existence  of  two  sequences  of  real 
eigenvalues  converging  to  ±°°. 

In  order  to  proceed  with  the  identification  of 
possible  accumulation  points,  we  make  the  important 
Assumption:   All  equilibrium  profiles  are  analytic 
functions  of  r ,  in  a  neighborhood  of  the  real  interval 
[0,R] . 

From  the  preceding  discussion  we  know  that  complex 

eigenvalues  can  accumulate  only  on  the  Alfven  or  cusp 

continua.   We  now  try  to  locate  the  accumulation  points 

in  the  continuous  spectrum.   Let  co  be  in  the  Alfven 

continuum,  and  let  r   be  such  that  A(r  ,co  )  =  0.   Let 

o  o   o 

also  A'  (r  ,  w  )  =}=  0  (where  the  prime  denotes  partial 
derivative  with  respect  to  r) .   From  the  implicit 
function  theorem  it  follows  that  there  exists  an 
analytic  function  r(co),  defined  in  a  neighborhood  of 
a)  ,  such  that  A(r(co),to)  =  0  and  r(u)Q)  =  rQ.   Moreover, 
r(w)  is  only  real  for  real  w(in  the  neighborhood  of  coQ) 
This  follows  from  the  identity 


|A   =  0  =  dr  A,  +  3A 
du>        dw       oo) 
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But  £ —  =  -2po  which  cannot  vanish,  otherwise  a  =  H  =  0 
3  01 

and  A'  =  0  contrary  to  our  assumption.   We  thus  have 

-t—  ±  0  ,  and  since  r  is  real  for  real  cu ,  it  cannot  be 

dco  ■ 

real  for  co  non  real. 

In  the  complex  r  plane,  equation  (B)  has  a  singu- 
larity at  r  =  r(co)  and  this  is  the  only  singularity  of 
the  equation  in  a  fixed,  small  enough  neighborhood  of 

r  ,  for  co  close  to  co  .   Let  co  have  a  positive   imaginary 
o  o 

part,  and  let  r(co)  have  a  negative   imaginary  part,  say. 

Let  £(r,w)  be  an  analytic 

solution  (in  co)  of  (B) 

which  satisfies  (B2) .  r-plane 

r 

The  value  £(R,co)  can 
be  obtained  by 


continuing  C  0  r  R 

o 


analytically  along 

a  path  r  in  the 

upper  half  r  plane  (see  fig.  1) 


*r(co) 
FIG.  1 


When  co-*co   equation  (B)  remains  regular  along  r,  hence 

£(R,co)  evaluated  this  way  is  analytic  in  co  and  cannot 

have  infinite  number  of  zeroes  near  co  .  A  similar  result 

o 

is  derived  for  co  in  the  lower  half  plane. 

This  process  shows  that  the  part  of  the  continuous 
spectrum  which  satisfies  the  assumption  above  (of  being 
a  simple  zero) ,  is  a  cut  in  the  domain  of  analyticity  of 

The  function  can  be  continued  analytically  through 
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the  cut  from  above  and  below.   A  similar  property  was 

derived  by  Tataronis  [8] . 

The  previous  argument  yields  similar  results  for 

a)   in  the  cusp  continuum.   We  are  left  with  only  two 

types  of  points  as  candidates  for  accumulation  points: 

i)    A  "tip"  of  the  continuous  spectrum,  which  is 

defined  for  the  Alfven  part  as  a  solution  oj  (with 

^  o 

corresponding  r  )  of  the  system  of  equations  A(r,oj)  = 

=  A'(r,oj)  =  0,  and  a  similar  definition  for  the  cusp 

continuum. 

ii)   The  point  of  overlapping  of  the  two  continua. 

At  this  point  A  =  C  =  0,  hence  H  =  a  =  0.   We  can  thus 

identify  this  point  as  follows:   Let  H (r  )  =0  (if  such 

a  point  exists),  then  to  =  -k*u(r  ).   In  this  treatment 

o    —  —  o 

we  exclude  the  possibility  H'(r  )  =  0. 

The  second  case,  although  included  in  the  first  one, 
is  singled  out  because  of  the  different  pattern  of  singu- 
larities of  the  differential  equation  near  r  .   In  case 
(i)  there  are  two  singularities  near  r   in  the  complex 
r  plane,  while  in  case  (ii)  the  equation  has  four 
singularities  originating  in  pairs  from  the  two  continua. 

Case  (ii)  reduces  in  the  static  case  to  H(r  )  =0, 

o 

oj   =  0  which  is  known  to  be  the  only  accumulation  point 
of  complex  eigenvalues.   It  will  be  shown  in  Chapter  IV 
that  a  similar  situation  occurs  in  the  steady  case.   The 
following  chapter  is  devoted  to  the  development  of  the 
mathematical  tools  required  for  the  analysis  of  our  case. 
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CHAPTER  III:   An  Asymptotic  Expansion 

III.l   The  Overlapping  of  Continuous  Spectra 

We  concluded  chapter  II  by  establishing  that  the 
point  of  overlap  of  the  Alfven  and  cusp  continua  is  the 
one  most  likely  to  serve  as  an  accumulation  point  of 
complex  eigenvalues.   The  asymptotic  method  to  be 
developed  in  this  chapter  is  tailored  for  the  investiga- 
tion of  this  case.   Although  it  is  possible  to  develop 
the  method  after  stating  general  assumptions,  we  prefer 
to  use  our  special  problem  as  an  explanation  for 
every  property  we  make  use  of. 

We  make  the  change  of  variables 


x  =  r-r    ,  c    =   w-u) 
o  o 


where  r   and  uj   correspond  to  the  overlapping  point  of 
o      o 

the  continua.   Without  loss  of  generality  we  assume  that 
x  ranges  in  the  interval  [-1,1].   We  rewrite  (B)  as 

(C)       [p(x,e)u']'  +  q(x,e)u  =  0   ,   -l£x<l 

The  proper  boundary  condition  at  x  =  -1  is  a  regularity 
condition  like(B2).   However,  all  we  use  is  the  analyti- 
city  of  the  solution  with  respect  to  the  parameter.   We 
therefore  can  assume  for      simplicity,  that  x  =  -1  is  a 
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regular  point  of  (C)  and  supplement  the  equation  with 
the  two  boundary  conditions 


(CI)  u(-l  ,e)  =  0 

(C2)  u(l  ,e)  =  0 


rAC 
we  recall  the  leading  coefficients  of  (B)  p  =   „  «■ 

q   a  +CKZ 

For  oj  =  tj  (e  =  0)  p  vanishes  with  its  first  derivative 

2 
at  r  =  r  (x  =  0) ,  thus  p  =  0(x  )  near  the  origin.   It 

is  also  seen   that  if  we  let 


(7)  z  =  - 


then  — »  p(ez,e)  converges  as  e-*-0  and  z  held  fixed   to  a 

e 
function  of  z  with  a  second  order  pole  at  infinity.   Like- 
wise we  find  that  for  e  =  0,  q(x,o)  approaches  a  constant 
as  x+0 .   Also  q(ez,e)  tends  to  the  same  constant  when 
first  e-*0  and  then  z->°°. 

Employing  these  properties  of  p,q  we  use  (7)  to 
rewrite  (C)  as  an  equation  in  z 

(D)  (p(z,e)u')  +  q(z,e)u  =  0 

P(z,£)  =  —  p(ez,  e) 
e 

q(z,  e)  =  q(ez,e) 
where  the  prime  denotes  differentiation  with  respect  to  z. 
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To  find  the  "eigenvalues"  e  of  (C )  we  use  a  boundary 
layer      analysis  similar  to  the  one  used  by  Pao  [9]  in 
his  investigation  of  the  static  screw-pinch  equilibrium, 
but  with  a  substantial  modification  required  by  the  pos- 
sibility that  the  eigenvalues  are  not  resticted  to  the 
imaginary  axis.   In  most  cases  asymptotic  methods  are 
believed  to  yield  correct  results  but  proofs  are  not 
available.   This  is  not  the  case  with  the  process  we  are 
about  to  describe.   The  proof  is  omitted  from  this  chap- 
ter for  the  sake  of  clarity,  but  will  be  presented  in 
much  detail  in  Appendix  A. 

We  divide  the  interval  [-1,1]  into  three  regions: 
(i)      The  outer  region  on  the  right,  in  which  x  Hei- 
di)    The  outer  region  on  the  left,  x  <  ~  I  £  I  • 
(iii)    The  inner  region  |x|<|e|' 

In  the  outer  regions  e  is  negligible  with  respect  to  x, 
while  in  the  inner  region  we  neglect  £  with  respect  to 
z.   To  leading  order  in  c,  an  eigenf unction  of  (c  )  is 
taken  to  satisfy  (c)  and  (C1-C2)  in  the  outer  regions, 
when  £  is  set  to  zero,  and  to  satisfy  equation  (D)  in 
the  inner  region  again  with  £=0.   The  matching  between 
the  different  regions  is  done  in  the  following  way. 
We  note  that  after  setting  c=0  in  ( C)  and  (D)  the  two 
equations  have  singularities  of  the  first  kind  (regular 
singularities)  at  x=0  and  z=^  respectively.   Moreover, 
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the  indicial  equations  that  determine  the  behavior  of 
solutions  near  the  regular  singularity,  are  identical 
in  both  cases  a   They  are  of  the  form 

(8)  v(v+l)  +  a  =  0 

where  the  real  constant  a  is  easily  derived  from  the 

coefficients  of  (C)  after  setting  e  =  0.   (8)  has  the 

two  solutions  v  and  -(v+1),  with  v  chosen  so  that 

Re  v  >  -  r  and  Im  v  >  0  when  Re  v  =  -  r  .   We  now 

require  that  the  outer  solution  with  behavior  x 

near  x  =  0 ,  be  matched  with  e~  times  the  inner 

solution  with  behavior  z       near  z  =  °°,  and  that 

xv  be  matched  with  e  [zv+cz      ]  where  c  is  some 

constant  determined  in  Appendix  A.   c  is  zero  in  the 

important  case  of  Re  v  =  -  ~     (v  complex) .   In  this 

last  case,  the  amplitude  of  the  eigenf unction  in  the 

1 
inner  region  is  e   2   times  larger  than  the  amplitude 

in  the  outer  region,  which  explains  the  name  "local 

mode" . 
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III. 2   The  Inner  Region  Equation 

The  coefficients  of  (D)  (after  setting  e=0)  can 
be  derived  from  (C)  by  substituting 


(9) 


H  -*■  hz 


o  ■*■   az  +  1 


where  h  =  H'(x=0)   ,  a  =  (k-u)'(x=0) 

All  other  functions  are  taken  to  equal  their  constant 

value  at  x  =  0 .   We  find 


(10)  p(z)  -  ^  A 

k^ 

2 
q(z)     =    -i^-    (BeH-pov)2    +    ^[pv2A-2paBQ     (Hv-B0a)]2 
rK   A 


2,1      o '  ,  '         2m, „.  . 

"r     (pP»    )       "    — 2-(Beh-pav) 

rK 


where  H  and  a   are  now  defined  by  (9) ,  r  =  r   and  as 


before 


0        7  2    2    2 

A  =  H  -po     ,    C  -  B*(BH  -p0  )  . 


The  inner  region  equation  now  has  coefficients  which  are 
rational  functions  of  z,  and  its  finite  points  of  singu- 
larity are  determined  by  the  vanishing  ol  There 
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are  four  such  points 


(11) 


'1,2 


h 

-a±  ,- 


'3,4 


/P 


-a±/6_  h 
/p 


The  point  z  =  °°  is  another  regular  singularity  of  the 
equation.  We  can  choose  two  independent  solutions  of 
the  equation  such  that 


(12) 


v±   =  zvf1(z) 


v2  =  .-«^«f2(« 


where  v  is  given  by  (8)  and  f  ,  f„  are  analytic  functions 
of  z  in  a  full  neighborhood  of  infinity,  and  satisfying 
f  (°°)  =  f_(°°)  =  1.   As  seen  from  (12),  the  two  solutions 
are  multivalued  functions  and  we  must  introduce  a  cut  in 
the  z  plane  that  extends  to  infinity.   We  choose  to  cut 
the  plane  along  the  negative  real  axis,  and  to  use 
-tt<  arg  z  <tt  in  (12)  . 


The  transformation 
(7)  maps  the  x  interval 
[-1,1]  onto  a  line 
segment  passing  through 
the  origin,  in  the  com- 
plex z  planes,  (see  fig.  2) 
Every  solution  of  (D)  along 


1 


z-plane 


FIG.  2 
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the  line,  can  be  continued  analytically  to  the  complex 
plane,  but  cuts  have  to  be  introduced  because  of  the 
presence  of  finite  singularities.   We  cut  the  plane 
along  the  real  axis  from  each  singularity  to  infinity, 
leaving  a  neighborhood  of  the  origin  in  the  domain  of 
regularity  of  the  solution. 

We  restrict  consideration  to  e  in  the  upper  half 
plane,  so  that 

(13)  0  <  0  =  arg  e  <  tt  . 

Transformation  (7)  then  maps  x>0  into  the  lower  half 
z  plane,  while  x  <  0  is  mapped  into  the  upper  plane. 
The  matching  of  innner  region  solutions  with  outer 
region  solutions  on  the  right  requires,  that  we  know 
their  behavior  near  infinity  in  the  lower  z  plane, 
while  the  matching  on  the  left  depends  on  the  behavior 
of  solutions  of  (D)  in  the  upper  half  plane.   In  order 
to  construct  a  solution  of  (C)  in  the  full  interval 
[-1,1] ,  we  need  to  know  how  solutions  of  (D)  change  their 
behavior  when  continued  from  the  upper  half  plane  to 
the  lower  one  in  the  cut  plane  of  figure  2.   It  is 
sufficient  to  know  how  v  ,v   of  (12)  change  their 
behavior.   We  describe  it  symbolically  as 

v      v    a    -(v+1) 

z   ►     +  Sz 

(14) 

z-(,+l)        .  5z-(v+l) 
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where  a,    S,  y,&   are  the  "connection  coefficients"  (not 
to  be  confused  with  the  same  symbols  used  for  different 
purposes  elsewhere  in  the  work) ,  and  the  arrow  points 
to  analytic  continuation  from  the  upper  to  the  lower 
plane.   z  ,  z       stand   for  v, ,  v  . 


The  connection  relations 
(14)  do  not  change  if  v..  ,  v 
are  continued  in  the  positive 
direction  around  all  the 
singularities  of  (D)  which 
are  in  the  right  half 
plane  (see  fig.  3).   This 
is  because  the  cut  in  the 


z-plane 


FIG.  3 


domain  of  regularity  of  the  two  functions  is  along  the 
negative  real  line  and  not  the  positive  one.   Clearly, 
if  (D)  has  no  singularities  in  the  right  half  plane,  then 
a=6=l,  B=y=0. 
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III. 3   Matching  the  Different  Regions 

We  begin  with  a  solution  u  of  (C)  (with  e  =  0)  in 
the  outer  region  on  the  right,  which  satisfies  (C2) . 
Near  the  origin  u  -  A(xv+rx"(v+  ')    with  r  a  constant 
and  independent  on  e .   xV,  x       stand  for  the  two 
solutions  of  (C)  with  the  corresponding  behavior.   u 
is  matched  with  a  solution  v  of  (D)  with  the  behavior 
v  -  A[cV(zV+cz"(v+1))+  rc-(v+1)z_(v+1)]  in  the  lower 
half  plane.   v  is  continued  into  the  upper  plane  and 
has  to  be  matched  with  a  solution  w  of  (C)  which 
satisfies  (CI)  and  thus  behaves  like  w  ~  a  (xV  +  £x       ) 
for  x  <  0  (arg  x  =  n)  ,  with  the  constant  I    independent 
on  e.   Using  (14),  matching  is  possible  if  and  only  if 

v     ,   v.-  -(v+1)    A   v 

ae      +  y(ce  +le  '    =   -  e 

a 


(15) 


v        v  „  -(v+1),    A,   v      -(v+1). 
Be   +  6(cc  +JU      )  =  -(ce   +  re      ) 


L<  ■  t 


2 

(16)  t  =  e 


then  (15)  can  be  written  as 


<17»  ♦  r(c  *  |) 


Bt  +  y(ct+  I)    =  -(ci  +r) 
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By  its  definition  Re  v  >  -  r  .   If  Re  v  >  -  —  then 
x  ■+  0  as  e  ■+  0 .   System  (17)  does  not  have  a  sequence 
of  solutions  t   converging  to  zero.   We  conclude  that  a 
sequence  of  complex  eigenvalues  e   -*-  0  can  exist  only  if 
Re  v  =  -  r  .   From  (8)  we  have 


v  =  J[-l+/T-4a] 


and  a  necessary  condition  for  the  existence  of  a  sequence 
of  complex  eigenvalues  near  e  =  0  is 


(18)  a  >  i 


In  the  static  equilibrium  case,  condition  (18)  is 
the  Suydam  condition  [10]  and  is  known  to  be  a  necessary 
and  sufficient  condition  for  local  instability  [11] . 

We  now  assume  that  relation  (18)  holds,  and  define 


(19)  u    =    /4a=T 


(20)  9  =  arg  e 


then  t  =  ellJ  =  exp(-u0  +  iy  log  |  e  |  )  .   Every  solution  t  of 
(17)  corresponds  to  a  sequence  of  numbers  e 


(21)      e   =  e      exp  (  -  — ^-)   ,   arg  e   =  arg  e   =6 
n    o^      v  ^      n     ^   o 


n  =  0,±1,±2, 
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with  6  determined  by 


(22)  |x|  =  e"M6 


As  n  ■+  °°  ,  e   ■+  0  which  is  consistent  with  the  assumption 
n 

underlying  the  asymptotic  expansion.   From  (21)  we  see 
that  asymptotically,  the  sequence  of  eigenvalues  is  on 

a  ray  6  =  constant,  with  a  spacing  of  exp ( )  between 

consecutive  members.   In  the  static  equilibrium  case  it 

2  TT 

is  known  that   e   is  real,  thus  6  =  s- 

Equation  (17)  is  equivalent  to  the  quadratic 
equation 


(23)         Bx2+U6  -  ra)x  -  r£Y  =  0 


where  we  used  the  fact  that  c  =  0  in  (17)  when  Rev  =    -   j   • 
It  is  seen  that  in  principle  two  sequences  of  complex 
eigenvalues  are  possible. 

The  following  property  of  r  and  His  sti  a, 

but  of  great  importance.   We  claim 

(24)  |r|  =  1   ,   |l|  =  e'^ 

proof:   Equation  (C)  has  real  coefficients,  and  so  arc  the 
real  boundary  conditions.  Hence,  in  the  ou-      ^jion  on 
the  right,  the  eigenf unction  is  a  multiple  of  the  real 
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solution  with  behavior 
x   +  rx 


2  ,  iy      ~iy. 
=  x   (x    +  rx   rj 


Since  we  chose  arg  x  =  0  for  x  >  0,  the  reality  of  this 
solution  requires  that  |r|  =  1.   The  proof  for  i    is 
identical  except  for  the  change  in  the  argument  of  x; 
arg  x  =  it  for  x  <  0 . 

In  view  of  (24)  we  define 


(25) 


ip       0     -yT7+iX 
r  =  e     ,    2,   =  e 


with  p,A  real.   We  now  prove  a  few  helpful  properties  of 
the  connection  coefficients. 
Property  1: 


det 


a 


=  1 


Proof:   Every  two  solutions  £-,r£2  of  (D*  (wnich  is  of 
a  self  adjoint  form)  satisfy  p(z)[£  ,£2]  =  constant, 
where  [^,,?2]  is  the  Wronskian  of  the  two  solutions. 
As  we  continue  the  solutions  around  the  singularities 
in  the  right  half  plane,   we  get  for  the  functions  (12) 


p(z)  [v1,v2]  =  p[aVl+6v2,  yv^&v^    =    (a6-3y)  Ptv!'v2] 
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Property  2 :   All  connection  matrices  associated  with  the 
same  path  of  continuation,  are  similar. 

proof:   Let  u  ,  v  denote  two  pairs  of  independent  solutions 
of  (D)  with  the  connection  relations 

u  •*■   Au   ,   v  ■*■   Bv  . 

We  can  express  the  solutions  u  as  a  linear  combination  of 
v:   u  =  Sv  with  a  constant  matrix  S. 

u  =  Sv  +  SBv  =  Au 


hence  B  =  S   AS 


property  3 :  a  and  6  are  real,y  =  -g 

Proof:   Equation  (D)  has  real  analytic  coefficients.   By  the 

Schwartz  reflection  principle  we  have 


p(z)  =  p(z)   ,    q(z)  =  q(z) 


Hence,  every  solution  v(z)  gives  rise  to  another  solution 


v(z)  which  is  its  reflection  in  the  real  axis,  v(z)  =  v(z) 
For  z  in  the  upper  half  plane  we  find 

v1(z)  -  aj  +  Hz  ' 
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thus  v..  =  $v.+av    Similarly  for  z  in  the  lower  half 
plane  we  get  v.  =  -yv.  +  ctv„ .   Comparison  of  the  two 
results  yields  the  desired  relation.   The  proof  for  6 
is  similar  when  we  consider  v„ .   It  is  noted  that  in 
the  proof  we  used  the  range   -tt  <  arg  z  <_  it. 

property  4 :   If  in  equation  (D) ,  the  coefficients  are 

2 
analytic  functions  of  z   rather  than  of  z,  we  have 

a+    6e    =  0 . 

Proof:   The  coefficients  of  (D)  are  real  on  the  imaginary 

axis.   Thus,  to  every  solution  v  corresponds  another 

solution  v  which  is  the  reflection  of  the  former  in  the 


imaginary  axis,  i.e.  v(z)  =  v(-z).   We  carry  now  a 

procedure  similar  to  the  one  used  before.   We  note 

—     ttI  — 
that  for  z  in  the  upper  plane  -z  =  e    z  ,  while  for  z 

in  the  lower  plane  -z  =  e    z.   The  difference  between 

the  two  cases  is  determined  by  the  choice  of  arg  z. 

Thus,  we  arrive  at  the  relations 


6(a+6e-M7T)  =  0    ,  |g|2  -  aV*  =  1 


from  which  it  follows  that  B  =f  0  and  a  +  6e~y7T  =  0 

From  equation  (23) ,  the  two  solutions  of  the 
quadratic  equation  satisfy  the  relation 


T1T2  =  "r£ 
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In  view  of  (24)  and  property  (3),  we  have  |t.t_|  =  e 
Let  |t,  =  e  1  (i=l,2)  ,  then  &2  =  it- 6,  and  we  have 
shown  the 

Result:   If  e=0  is  an  accumulation  point  of  complex 
eigenvalues  ,  then  they 


-UTT 


e-plane 


converge  (asymptotically) 
along  two  rays,  which  are 
the  reflection  of  each 

other  in  the  imaginary  .. 

/  I 

axis  (see  figure  4)  .  * 

* 
(The  rays,  not  the 

eigenvalues,  are 

reflected) .   The  FIG-  4 

eigenvalues  in  the  lower  half  plane  are  the  complex 

conjugate  of  those  in  the  upper  half  plane. 

It  was  assumed  all  along  that  0<6<tt.   As  was  pointed 

out  before, this  range  of  arg  e  was  used  in  determining  the 

connection  coefficients  as  well  as  in  deriving  some 

properties  of  them.   Thus,  only  solutions  of  (23)  which 

correspond  to  this  range  of  0,  give  rise  to  local 

instability.   In  order  to  find  when  this  happens  we 

rewrite  (23)  as 


B 
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Using  (25)  and  defining 
(26)  t  =  s  exp  y(P+*-2b) 


B  =  |6|eib   ,   *  =  2ZL 


we  have 


1    .    i*       -MTT    -ilj,. 


(27)     s  +  - -  -i-  UeXT  -  e    6e  xr] 


The  angle  of  convergence  9  can  be  determined  by 
si  =  It  I  =  e"M6.   Let 


-yrr 
F(s)  =  s  + 


(28) 


^  =  TrT    ~  e 

The  mapping  F(s)  transforms  the  circles  in  the  s  plane 
that  are  centered  at  the  origin,  into  a  family  of  confocal 
ellipses  in  the  F  plane  with  foci  at  ±2  exp(-klJL).   In 
particular,  the  circle  |s|  =  e    '  *"   maps  onto  the  line 
segment  (degenerate  ellipse)  between  the  two  foci.   The 
transformation  G(i|>)  maps  the  real  line  onto  a  horizontal 
or  vertical  ellipse  (to  be  denoted  by  a)  with  major  axes 
1     a  +  e_1JTr6  I  "  This  ellipse  is  not  a  member  of  the 
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former  confocal  family.   Indeed,  belonging  to  this 

2 

family  requires  that  a6+ | 6 |   =  0,  a  contradiction. 

We  arrive  at  the 

Result:   As  ty   varies,  the  solutions  of  (27)  change 
their  absolute  value.   Thus  the  angle  of  accumulation 
of  complex  eigenvalues  depends  on  the  global  behavior 
of  the  equilibrium  profiles  and  not  only  on  their  value 
at  x  =  0 .   In  particular,  the  angle  will  change  if  we 
change  the  position  of  the  outer  wall. 

This  situation  differs  from  the  static  equilibrium 
case,  where  the  self  adjointness  of  the  problem  provides 

2  IT 

that  e   is  real  and  9=2   ■   A  derivation  of  this 
result  based  on  our  formulation  is  obtained  as  follows. 
When  the  flow  vanishes,  the  coefficients  of  (D)  depend 
on  z   only.   From  property  4  we  get  a+6e    =  0  ,  and 
the  a-ellipse  degenerates  to  a  horizontal  line  segment. 
Equation  (27)  then  yields  |s|  =  e  £■   regardless  of  ty . 

Proceeding  with  the  investigation  of  (27),  we 
require  the  solutions  to  satisfy  e~yTT<|s|<l,  corresponding 
to  tt>0>O.   The  point  G(ij/)  must  then  be  in  the  interior  of 
the  ellipse  corresponding  to  |s|  =  1  (and  to  | s |  =  e    ), 
which  is  characterized  by  its  major  axes  1  i  e 
We  denote  this  ellipse  by  u.      Three  cases  seem  possible 
(see  figure  5) . 
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(a) 


(b) 


FIG.  5 
case  (a):   a  is  in  the  interior  of  y,  which  yields 

eigenvalues  for  all  \\> . 
case  (b) :   a  intersects  y.   There  exist   eigenvalues 

for  some  ^  but  not  for  all  values, 
case  (c) :   a  is  in  the  exterior  of  y.   There  are  no 

eigenvalues  near  e  =  0. 
In  view  of  properties  1  and  3,  the  connection 
coefficients  satisfy  the  relation 
(29)   a  6  +  |3|2  =  l 

which  we  use  to  express  | 3 |  in  terms  of  a  and  6 .  The 
three  cases  can  now  be  characterized  in  the  following 
way.   Let 


(30)     Q(ct,6)  =  a2eU7T  +  2  cosh  (yn)a6  +  62e~M1T 


2  yrr 


case  (a)  corresponds  to  Q(a,6)<  4  sinh  *—*■ 
(31)   case  (b)  :   4  sinh2  ]p<  Q(ct,6)<4  cosh2  ^L 


case  (c) 


Q(a,6)  >  4  cosh 


2  yjT 


» — •  ' 
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The  derivation  of  these  inequalities  simply  uses 
comparison  of  the  major  axes  of  the  two  ellipses  a,u. 

A  further  simplification  in  the  characterization  of 
the  three  cases  (31)  is  required  in  order  to  enable 
us  to  determine  which  one  of  them  holds,  while  bearing 
in  mind   that  the  connection  coefficients  a, 6  are  not 

known  explicitly.   We  recall 

z-plane 
property  2  which  states  that 

r2  ri 

all  connection  matrices 

associated  with  a  path 

of  analytic  continuation 

are  similar.   In  particular, 

FIG.  6 

their  trace  is  independent 

of  the  pair  of  independent  solutions  we  continue. 

Let  tr,  ,  tr   be  the  traces  of  the  connection  matrices 

associated  with  the  paths  r\,  T_  winding  around  the 

singularities  of  (D)  in  the  right  and  left  half  planes 

(see  figure  6).   We  use  the  two  solutions  (12)  to  obtain 

the  following  relations 


(32)  a+6   =  trl 


aeM   +  6e     =  -  tr2 
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from  which  a,  6  can  be  determined  in  terms  of  tr..,  tr_ 
The  factors  in  the  second  equation  of  (32)  arise  from 
the  dependence  on  arg  z  in  the  definition  of  (12) , 
which  changes  the  connection  relation  of  v..  say,  to 


v     ,   -2iTi.v    Q/   -2iri  -(v+1) 
z   ■*■  a  (ze     )   +6  (ze     ) 


when  v,  is  continued  along  T2    in  the  negative  direction. 
The  direction  of  continuation  however,  is  immaterial  for 
our  purposes.   The  traces  of  a  matrix  and  its  inverse 
are  equal  if  the  determinant  is  unity  (see  property  1) . 
We  observe  that  by  (32)  we  have 

(33)  Q(a,6)  =  -tr1-tr2   . 

If  in  particular  the  two  traces  (which  by  (32)  are  real) 
have  equal  signs,  then  only  case  (a)  of  (31)  holds  and 
eigenvalues  are  possible  regardless  of  the  global  behavior 
of  the  plasma.   We  will  pursue  this  point  in  chapter  IV 
and  show  that  this  is  indeed  the  situation  for  many 
plasma  profiles. 

Analytic  determination  of  the  two  traces  generally 
requires  that  we  know  explicitly  two  indpendent  solutions 
of  (D)  (e.g.  in  the  form  of  an  integral  representation). 
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Hov;ever,  if  one  of  the  half  planes  contains  only  one 
regular  singularity,  the  associated  trace  can  be 
determined  by  the  characteristic  exponents  of  solutions 
near  the  singularity.   For  the  indicial  exponents 
\.,    \2,    the  trace  is  e27T1  Al  +  e27T1  A2  .   If  a  numerical 
determination  of  the  connection  coefficients  is   sought, 
relation  (32)  offers  the  great  benefit  of  avoiding  the 
necessity  to  compute  to  infinity,  especially  when  the 
solutions  of  (D)  are  oscillatory. 

We  now  proceed  with  proving  the  asymptotic  nature 
of  the  matching  procedure  described  in  this  chapter, 
after  which  we  will  resume  the  investigation  of  the 
screw  pinch  steady  equilibrium. 
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Appendix  A 

The  asymptotic  expansion  used  in  chapter  III,  will 
be  carried  out  rigorously.   We  first  describe  the  behavior 
of  solutions  of  (c)  in  the  outer  and  inner  regions 
associated  with  the  singularity  at  x  =  0 ,  e  =  0. 
Equation  (c)  has  the  form 


(Al)      (p(x,£)u')'    +   q(x,e)u  =  0,   -1<x<1 


Let  u(x,e)  be  a  solution  of  (Al)  which  behaves 
like  x   near  the  origin,  when  e  =  0.   We  require  that 
u  satisfies  two  initial  conditions  at  x  =  -1  independently 
of  e,  such  that  it  has  the  desired  behavior.   Let 
<J)  (x)  =  u(x,0)    and  define  v(x,e)  such  that 


(A2)  u  =  <f>v 


v  then  satisfies  the  equation 


(A3)         (p<f>2v')'  =  [Aq4>2  +  <MAp<j>')  ']v 


with  initial  conditions 


(A4)  v(-l,e)  =  1   ,   v'  (-l,e)  =  0 
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where 


(A5)  Aq  =  qQ-q    ,    Ap  =  PQ  "  P 


with  qQ  =  q(x,e=0)   and   pQ  =  p(x,e=0) 


We  assume  that  <$>   does  not  have  any  zeroes  in  -1<*<0  , 

otherwise  we  impose  initial  conditions  for  u  and  <j>  at 

-l<x  <0.   The  right  hand  side  of  (A3)  can  be  viewed  as 
o 

an  inhomogeneous  term,  and  it  is  seen  that  the  homo- 
geneous equation  has  the  two  independent  solutions 


-=±?    .   v  therefore,  satisfies  the  integral 
J    p<T 


1. 

equation 


(A6)    v(x,c)    =    1+         K(x,S)  [Aq<p2+<p(Ap<p')  ']v(£)d£    =    1+Kv 

-1 

where 


(A7)  K(> 


dt 
P9     (x) 


r  AC 
We  recall  that  p  =  -=— r y   -   (See  the  definitions  in 

P/a*+CK 
section  1.2). 

For  e  and  x  small,  A  has  two  zeroes,  at  x.  ,  =  h.  , 
with  h.  (i  =  l,2)analytic  functions       ind  h^O)  =  0. 
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We  can  thus  write  p  as 


p  =  f  (x,e)  (x-h1(e))  (x-h2(e))  . 


We  restrict  our  treatment  to  |x|>N|e|  ,  |e|<6  with  N 

large  enough  such  that  the  four  singularities  near 

x  =  0  are  in  the  interior  of  |x|  =  H|e|  for  |e|<6. 

Despite  the  denominator  of  p,  f  (x,e)  can  be  uniformly 

bounded  0<c1<|f|<c2  for  N|e|<|x|<l  ,  |e|<6  . 

The  integral  equation  (A6)  is  solved  by  the  Picard 

iterations  method.   We  define  v    , ,  =  v   +Kv»   v   =1 

n+1     o      n  '    o 

The  integral  in  £  is  taken  along  a  path  composed  of 

Y-j-the  real  interval  -1<£<-|x|,  and  y2~  the  circular 

arc  |x|  =  const,  from  -|x|  to  x  in  the  negative 

direction  (see  fig.  Al) .   Proceeding  in  the  standard 

way,  we  look  for  a  bound  of  K 

which  will  guarantee  the  x-plane 

convergence  of  the 

iteration  scheme. 

In  this  appendix  we 

use  the  notation 


K(x,5)  =  K(z,0  = 
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where  <J>  =  x  g(x),  and  g  is  bounded  away  from  zero.   For 
| t | >N  ,  the  quadratic  in  the  denominator  is  close  to  t 
(if  N  is  so  chosen)  ,  and  we  get  (for  \>   real) 


|K(zf  C)  I  < 


const 

.  I  2v+l 


+       V*     i 


2v+2    |  ,2v+l 
z 


The  last  term  is  the  contribution  of  the  circular  arc 
We  have 


c         c 

(A8)  |K(z,c)|  <  5— y  [  2^TT  +  T^l   ]       • 

|e|2v+l    |z|2v+l     |C|2^+1 


From  their  definition 


Aq  =  0(±) 


Ap  =  0(ex) 


The  two  kernel  terms  <p    Aq  and  <J>(Ap<f>')  '    are  of  similar 
order  when  we  fix  z  and  let  c  be  small,  and  it  is 
sufficient  to  treat  one  of  them. 


K(x,OAq<{>  U)<H|.    const 


2v-l 


Yl 


!zl 


+    Be       +    c    £  £ 


r  _1  +  2        , 

|z|2v+l  |;|2V+1JUI  d|C 


.-(2v+l)    with    constants    A/B,C, 


50 


The  integral  along  Y~  has  a  bound  of  ,  and  can  be 

|z| 
incorporated  in  the  previous  result. 

x 

Let         F(x,e)  = 


|K(X,C)  [Aq4>2  +  ^Ap^')  'del 


-1 

Fn+1 

then  |v  ,,-v   <  -> — m-  ,  as  can  be  seen  by  induction.   We 
1  n+1   n1  —  (n+1) 

write  v   =  v  +v  ,-v  +...+V  -v   ,    .    The  sequence  converges 
n     o  l  o      n  n-1  * 

uniformly  in  a  compact  set  of  z.  Let  v  ■*■  v,  then  v 
satisfies  the  integral  equation  (A6),  hence  also  the 
differential  equation  (A3) ,  and  we  have 


(A9)     |v-l|  <  eF-l   ,   F  <  -A_  +  Bl  el  +  clel  ieZr(2v+1) 


if  -(2v+l)>0   ,  i.e. 


(A10)  v  <  -  j 


(which  we  now  assume,  although  it  is  contrary  to  previous 
notation)  then  for  fixed  z,  F  remains  bounded  as  e+0 . 

We  want  to  show  that  as  e-»-0  and  z  kept  fixed,  v(z,e) 
converges  to  a  function  v..  (z)  which  is  the  solution  of 
(A3)  with  e  =  0.   ((A3)  has  to  be  multiplied  by  e~  V 
in  order  for  the  coefficients  to  converge  to    non  trivial 

functions  of  z)  .  From  (A6)  we  get 


dv     1 


dx     ,2,  . 
P4>  (x) 


[Aq<}>2  +  4>(Ap<J>')  ']vU)d£ 


-1 
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Hence 

-1/e  A 


dv I     <      const 
STi"1    -    I     i2v+2 


?2Mv(Od?    I    <   const   -^[Iz^^cjer2"] 

I  z  I 


which  proves  that  -^ —  remains  bounded  as  e  ■*■  0  . 
c  dz 

By  the  Arzela-Ascoli  Theorem  there  exists  a  subsequence 
v(z,e  )  converging  uniformly  to  v, (z)  as  e_~*"0 •   Because 
of  the  uniform  convergence,  the  limit  function  v. 
satisfies  (A3)  when  z  held  fixed  and  e  =  0.   Equation 
(A3)  with  e  =  0  has  two  independent  solutions  which  behave 
like  1  and  z~      v+1>    near  z  =  °°  ,  of  which  the  latter  is 
the  larger  one  by  (A10) .   (A9)  shows  that  v   behaves  like 
1  near  infinity  and  thus  is  exactly  the  smaller  solution. 
Since  we  can  identify  the  limit  function,  it  follows  that 
v  (z  ,e)  -*-v,  (z)  for  any  sequence  of  e.   (Otherwise,  there 
exists  a  sequence  v(z,e.)  which  stays  away  from  v,(z), 
but  a  subsequence  of  it  must  converge  to  v.) . 

From  the  definition  u  =  <f>v  we  find  that  as  e-*0, 
e~Vu-»-zvv1  (z)  ~  zv  where  u  is  the  large  solution  near 
x  =  0(v  <_   -  j)     .       In  the  important  case  of  Re  v  =  -  y 
and  v  complex,  the  proof  follows  throught  for  the  two 
independent  solutions. 

The  small  solution  u  -  x       requires  a  similar 
treatment.   We  use  the  same  iterative  scheme,  but  this 
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time  with 


(A10) 


v   =  -  (2v+l) 
o 


^     x  c2v+1 


Ne 


P<J> 


with  a  path  of  integration  consisting  of  a  circular  arc 
and  a  ray,  (see  fig.  A2 ) .   For  e  =  0,e       vq   ~  x 

x-plane 
x 

near  the   origin  and  in  fact, 

it  converges  to  this  function 

as  e   ■*■   0  and  x  kept  fixed. 

(Dominated  Convergence 

Theorem) .   Carrying  out 

the  iteration  process  we  FIG.  A2 

find,  that  the  limit  function  v  exists  and  satisfies 


-1 


\ 


(All) 


Iv-v0|  1  (e*-l)  |vo| 


For  fixed  z 


=  -(2v+l) 


dx 


N 


fgx 


2v(t-  ^)(t-|2.) 


and  as  e  ->  0  it  converges  to  a  function  of  z  with  a 
behavior  of  z"(2v+1)  near  infinity.   Again,  using  the 
Arzela  Ascoli  theorem,  there  exists  a  sequence  en  ■+   0 
such  that  v(z,£  )  -*■  v2(z).   By  (All),  v2(z)  -  z 
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near  z  =  °°  and  therefore  v_  is  the  sum  of  this  solution 

of  (A3)  (with  e  =  0)  and  cv.  with  a  constant  c,  and  v.. 

as  defined  before.   The  constant  c  is  determined  by  the 

requirement  that  v_(z=N)  =  0,  which  follows  from  (All) 

and  v  (z=N)  =  0.   The  unique  identification  of  the 
o  ^ 

limit  function  enables  us  to  conclude  that  v  (z  ,  e) -+-V- (z) 
for  any  sequence  e    -*•  0 . 

It  is  important  to  note  that  the  iteration  procedures 
used  before,  also  converge  when  x  is  considered  the 
argument  and  e  is  small.   As  e  ■*■  0 ,  the  two  limit 
functions  converge  to  1  and  to  a(x) «x        with  a(0)  =  1, 
respectively. 

We  now  prove  the  validity  of  the  connection  process 
used  in  chapter  III.   We  define  u.  (x,e)  ,  u„(x,e)  in  the 
neighborhood  of  x  =  -1,  as  the  two  functions  obtained 
by  the  iteration  procedures,  multiplied  by  <?>(x).   These 
functions  are  continued  analytically  through  the  upper 
half  plane  above  the  singularities 
of  equation  ( C)  (see  fig.  A3). 

We  use  these  solutions  to  solve  x-olane 

the  initial  value  problems  at 
x  =  ^1  (in  the  outer  regions)  ,     , 


-1  0 

and  obtain  solutions  that  can  • 

j      ,       j.0       ,  FIG-  A3 

be  expressed  as  a(u_  +  f-u.)  , 

A(u.+ru  )  on  the  left  and  the  right,  with  a,  A,r,  I 


54 


depending  on  e  only.   It  is  noted  that  u,  (x,e)  is 
analytic  with  respect  to  e,  but  that  u~(x,e)  is  only 
continuous  in  c    (for  v  real)  because,  by  its  definition, 
v   in  (A10)  is  not  analytic  in  e.   However,  continuity 
when  Re  v  4  ~  t"  suffices  for  our  purposes.   When  u,  ,  u~ 
are  continued  along  the  real  axis  from  x  <  0  to  x  >  0 , 
we  express  the  resulting  functions  in  terms  of  the  values 
of  u- ,  U2  previously  obtained,  i.e.  by  connection 
coefficients.   These  coefficients  have  simple  dependence 
on  the  connection  coefficients  a,3»Y/<5  of  the  functions 
e   u,  ,  e    u2  *   T*ie  -I-ast  functions  converge  as  e  -*■  0 
and  z  is  held  fixed,  and  therefore  the  connection 
coefficients  associated  with  them,  converge  as  e  ■+■  0  to 
the  values  that  were  used  in  the  asymptotic  treatment  of 
chapter  III.   Finding  complex  eigenvalues  requires  that 
we  solve  the  equation 


(A12)  3t2  +  U6-ra)x  -  r£y  =  0, 


where  x  =  e    and  all  the  coefficients  are  analytic 
functions  of  e .   (It  is  clearly  necessary  that  v  be 
complex.)   If  instead  of  the  coefficients  in  (A12)  we 
take  their  value  at  e  =  0  (and  keep  the  same  notation) , 
equation  (A12)  can  be  written  as 


gT2  +  (£6-ra)x  -  rly   =   0(e) 
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By  a  Rouche  theorem  argument,  it  is  easy  to  show  that 
the  homogeneous  and  nonhomogeneous  equations  have  the 
same  number  of  zeroes . 
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CHAPTER  IV:   The  Generalized  Suydam  Condition 
IV. 1   The  Point  of  Overlap 

We  now  use  the  results  obtained  in  chapter  III  to 
decide  the  question  of  local  stability  near  the  point 
of  overlap  of  the  Alfven  and  cusp  continua.   The  inner 
region  equation  associated  with  this  point,  has  four 
finite  regular  singularities,  given  by 


(ii)     zlt2  =   -^r-   ,  z3 


-a±  -h       J'     -a±/g  h 


Assuming  for  simplicity    a  >  0,  three  situations  are 
possible : 

2    2 

i)       h  -pa   <  0  ;  All  four  singularities  are  in  the 


(34) 


left  half  plane. 


2      2     2 
n)      Bh   <  pa   <  h   ;   Only  z,  is  in  the  right  half 

plane. 

2      2 

iii)     3h   -  pa   >  0  ;   zi'zo  are  i-n   tne  right  half 

plane,  the  other  two  points 
are  in  the  left  plane. 

(a, 3  here  are  defined  in  (9)  and  are  not  connection 
coefficients) . 
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If  relation  (34. i)  holds,  it  was  indicated  before 
that  the  connection  matrix  is  the  identity  matrix,  and 
a  sequence  of  complex  eigenvalues  cannot  exist  near 

E  ■  0. 

Next  we  deal  with  case  (34.iii).   From  expression 
(10)  we  find  the  behavior  of  solutions  of  (D)  near  the 
singularities.   Near  the  Alfven  singularity  z,  (  and 
likewise   z  )  two  independent  solutions  can  be  chosen 
to  be 


(35) 


u,(z)   ,   u1(z1)  =  1 


u_(z)  =  u,(z)  log  (z-z1)+  f(z) 


with  u  ,f  analytic  functions  in  a  neighborhood  of  z  . 
Near  the  cusp  singularity  z  (and  likewise  z 4) ,  there 
exists  two  solutions  of  the  form 

(36)  v1{z)        ,   vx(z3)  =0   ,   v{ (z3)  =  1 


v2(z)  =  v^z)  log  (z-z3)  +  gvz) 


with  v  ,g  analytic  near  z^  . 
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Let  (ul)  =   S (  1)   with  S  =  (    ,)  a  constant  matrix. 
u2        v2  c   d 

We  now  continue  the  solutions  (35)  along  a  path  which 
winds  around  the  two  singularities  z  ,  z^  (see  fig.  7) 

Clearly  <«1,  ♦  ,^  J,  s  (^  J,  B^l)  . 

A  simple  calculation  shows  that 

the  trace  of  this  connection 

matrix  is  z -plane 

(37)   tr,  =  2  +  47T  *> 


1  A 

where  A  is  the 


Z4  Z2 


determinant  of  S  and 

b  is  its  upper  off  FIG.  7 

diagonal  element.   By  the  definition,  u.  =  av,  +  bv   . 

[u1,v1]  X  l  l 

Hence  b  =  -f ^-  with  the  brackets  again  standing  for 

lv2'vlJ 
the  Wronskian.   By  the  choice  of  the  solutions,  b  is 

2 
real  and  therefore  b   is  non  negative.   As  to  A, 

[u.  ,u_]  =  Atv-j/V,,].   Since  equation  (D)  is  of  a  self 

adjoint  form,  p[u.u„]  =  constant  for  any  two  solutions 

u..  ,  u_  .   We  write 


(38) 


p[u17u2] 

A   =  T  ~ 

p[v1,v2] 


and  evaluate  the  numerator  for  z  -»■  z,  ,  while  the  denomina- 
tor is  evaluated  at  z  ■*■   z  _  . 

2 
ul 
tul'u2]  =  [ul'f]  +  i= 


Zl 
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By  (10)  p  =  ^yA(z)  .   As  z  ■*  z  ,  p  ~  -jA'  (z,  )  (z-z,  ) 
K  K 


and  p[un,u  ]  =  —*■  A'  (z.)     .   A  is  quadratic  in  z  with 
1      Z  K       i 

2    2   2 
a  leading  term  (h  -pa  )z   ,  and  z,  is  its  rightmost 

zero.   Hence  A'(z  )  >  0  and  the  numerator  of  (38)  is 

positive . 

2 

vl 


[v1,v2]    =    [v1,g]    +   — 


Z3 


By  (36)  p[v,,v  ]  =  -pv  'g(zj.   g  satisfies  the  equation 


pv. 


i 


CpgT  +  qg  =  -2^-    -   (^f-J'v, 


and  is  a  regular  function  near  z^.   If  qQ  is  the  coefficient 
of  —2 —  in  the  expansion  of  q  ,  we  get  q  g(zj  =  -pv  '(z3) 

2 
or  -pv'g(z-)  =  q  g  (z,).   As  seen  from  (10),  the  sign  of 

q   =  sign  AC'(z^)  which  is  clearly  positive.   Therefore 

in  (38)  A  >  0  and  we  have  tr   >  0.   We  now  repeat  a  similar 

process  for  the  evaluation  of  tr  .   The  same  argument 

follows  through,  where  we  have  to  decide  the  sign  of 

A'  vz2)A(z4)C  (z4)  .   Clearly  A'(z2)  <  0,  A(z4)  >  0, 

c'(z  )  <  0  ,  and  again  we  get  tr2  >  0.   Case  (34.iii)  is 

therefore  in  the  category  of  (31. a),  and  accumulation 

of  complex  eigenvalues  is  possible  regardless  of  the 

position  of  the  boundary  or  the  global  properties  of  the 

plasma.   Case  (34.  ii)  is  not  as  simple  as  the  others  and 
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its  full  details  are  yet  to  be  worked  out.   Although 
in  the  right  half  plane  we  can  use  the  two  solutions 

(35)  to  determine  easily  that  tr   =  2/      in  the  left 
half  plane  we  have  three  singularities  and  the  previous 
methods  fail  to  yield  the  sign  of  tr„.   It  is  instruc- 
tive to  mention  that  at  least  in  some  special  cases, 

(34. ii)  does  yield  results  similar  to  (34.iii).   In 

the  limit  3  ■+  0  (zero  speed  of  sound)  the  two  cusp 

singularities  coalesce  into  a  regular  singularity 

located  at  z  =  —  in  the  left  hand  plane.  If  further- 
ex  ^ 

more  the  equilibrium  flow  has  no  azimuthal  component 
(v  =  0),  this  point  is  a  regular  point  of  the  equation. 
Accordingly,  the  only  singularity  in  the  left  half 
plane  is  z„  and  tr~  =  2.   Because  of  the  continuous 
dependence  of  tr   on  3  and  v,  tr„  will  remain  positive 
for  3  and  v  small  enough,  so  that  (31.a)  holds. 

It  would  be  useful  to  work  out  a  special  case  in 
detail,  and  to  convince  ourselves  of  the  added  feature 
of  non  self  adjoint  systems,  namely  the  global  depen- 
dence of  the  angle  of  accumulation.   The  case  we  describe 
is  incompressible  flow,  that  is,  a  divergence  free 
perturbed    velocity  field.   This  case  can  be  reduced 
formally  from  the  Lundquist  equations  (L)  by  letting 

the  speed  of  sound  approach  infinity,  which  amounts  to 

2 
letting  B*  -*■  °°    ,  p  -*■   1.   Proceeding  with  the  asymptotic 
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expansion  described  in  chapter  III,  we  find  that  the 
inner  region  equation  is: 


(39)  (Au1)  '  +[|  ^(BQH-pav)2-rK2(p  p°  '  )  '-  f^d^n-P  va  )  ]  u  =  0 

r 


The  term  containing  the  two  cusp  singularities  disappeared, 

and  we  have  only  three  regular  singularities  in  the 

plane:   the  point  at  infinity  and  the  two  Alfven  singulari- 


ties  z 


1,2      M/ r 
-a±  //£ 


which  are  in  different  sides  of 


2    2 
the  imaginary  axis  if  h  -pa   >  0  .   Assuming  this  inequality 

to  hold,  the  indicial  equation  that  determines  the 


behav 


ior  (zc,  say)  of  solutions  near  z,  2  i; 


If  we  let 


2     \       4k  2 

A  2   r2   6 


1,2 


1,2 


|^-(BeH-pav) 


=  0 


k(B0+/p~  v) 


;1,2     r(h+/p  a) 


x  =  0 


then  tr.  2  =  2  cosh(2TT  X^^  2)  . 

By  (32)  we  find  that  a  +  6e~y7T  =    —  sinhn  ( X, +X0 )  sinhTT  (  X, -X0) 


1+e 


pTT 


.    - 


;  2 


This  expression  vanishes  (corresponding  to  a  90°  angle  of 
accumulation)  if  and  only  if  1\     ,    which  clearly 

is  not  the  general  case.   It  should  be  noted  that  the 
incompressible  limit  is  a  regula  r  limit  for  the  connection 
coefficients,  and  that  as  a  result  we  can  reach  a  similar 
conclusion  for  a  mgh.   The  reason     hat  the 
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connection  matrix  is  derived  from  solutions  of  the 

equation  along  a  finite  path  which  is  bounded  away 

2 
from  any  singularity  as  a   -*■<». 
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IV. 2   The  Tip  of  the  Continuum 

It  was  shown  at  the  end  of  chapter  II  that  beside 
the  point  of  overlap,  a  tip  of  the  continuum  is  the  only 
possible  accumulation  point  of  unstable  modes.   We  first 
consider  a  tip  of  the  Alfven  continuum,  given  by 


A(r  ,u>  )  =  A'  (r  ,u  )  =  0 
o  o        o  o 


Let  x   =  r-r   ,  then  near  x  =  0  and  ifor  cj  close  to  w 
o  o 


2    2     2    2  2     2 

A  =  H  -pa   =  H  -pa  -2pa  (w-io  )+0(u)-u)  )  -  ax  -2po  (oj-oo  ) 

o  K  o    o       o  o    o 


where 


o   =  o  (x ,  a)=u)  ) 
o  o 


a  =  j   A"  (x=0,  uj=io  ) 


If 


2po 


>  0  ,  the  part  of  the  Alfven  continuum  with  the 
x=0 


tip  oj   is  in  the  half  line  u>  >  u>  .   For  simplicity  we 

c      o  —  o 

assume  this  to  be  the  case  and  we  define  e  such  that 

e   =  uj-uj   and  Im  e  >  0 .   If  —£2.  <  0  we  set  c   =  oj  -w. 
o  —         a  o 

Proceeding  as  in  chapter  III,  we  divide  the  interval  of 
x  into  outer  and  inner  regions  corresponding  to  |x|  >  |e| 
and  |x|  <  |e|.   Using  equation  (B)  to  describe  the  plasma, 
we  find  that  det  A  vanishes  when  A  vanishes,  and  therefore 
when  x  and  e  are  small  enough,  the  outer  region  i      on 
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takes  the  form 


(40)  [(x2  -  ^££  e2)u']  '  +  qu  =  0 


with  the  coefficients  evaluated  at  x  =  0 ,  e    =  0 .   In 

the  inner  region  we  define  z  =  —    ,  and  for  z  fixed 

e 

and  e  -*■  0 ,  we  get 

,    2pa 
(41)  [(zZ  -  — )   u'] '  +  qu  =  0 

a. 

which  is  Legendre's  equation.   In  contrast  with  the 
overlapping  case,  this  equation  has  only  3  singularities 
in  the  complex  z  plane  and  not  five  of  them.   As  often 
happens  with  special  functions,  we  do  not  need  the  avail- 
able explicit  solution  of  equation  (41) ,  but  derive  the 
required  results  from  the  equation  itself.   As  in  (18)  , 
a  necessary  condition  to  have  an  accumulating  sequence 
of  complex  eigenvalues  e  near  c   =   0 ,  is  that  q  >  -r   . 
Introducing  the  connection  coefficients  as  before,  we 

find  that  ct+6e_lJ71  =  0,  since  the  coefficients  of  (41) 

2 
depend  on  z   alone.   We  thus  find  that  asymptotically,  the 

complex  e  eigenvalues  are  purely  imaginary,  corresponding 

to  modes  oj  <  u   on  the  real  axis.   As  assumed  before,  this 
o 

ray  is  not  part  of  the  continuum  and  we  have  established 
that  eigenvalues  converging  to  u  must  be  real.   This  result 
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is  not  only  asymptotic  but  exact,  as  can  be  derived 

from  equation  (B)  by  an  oscillation  argument  as 

follows.   For  oo  =  co   the  equation  is  singular  at  x  =  0 . 

Also,  if  q  >  —   the  solutions  oscillate  rapidly  near 

the  origin.   Thus,  for  co  close  to  oo   and  real,  but  not 
3  '  o 

in  the  spectrum  (oo<oo  )  ,  every  solution  of  (B)  has  many 
zeroes (although  not  an  infinite  number)   near  x  =  0, 
which  is  a  consequence  of  the  analytic  dependence  of  the 
equation  and  its  solutions  on  oj  ,  when  we  consider  |x|>6>0 
with  6  small  as  we  wish.   More  precisely,  we  consider 
the  solution  that  satisfies  u(-l)  =  0  ,  u'(-l)  =  1. 
As  we  let  co  -»■  oo   ,  zeroes  of  u  appear  at  x  =  1  and  move 
into  the  interval  (-1,1).   Every  such  zero  yields  an 
eigenvalue  and  the  existence  of  many  zeroes  of  the  so- 
lution in  the  range  of  the  plasma  for  co  close  to  oo  , 

o 

indicates  the  existence  of  a  sequence  of  eigenmodes 

ca>   ■*  ca)   .   The  existence  of  a  sequence  of  real  eigenvalues 
no 

converging  at  a  tip  of  the  continuum  was  discovered  by 

Tataronis  in  the  static  equilibrium  case.   It  was  then 

shown  [12]  to  be  either  Sturmian  or  anti  Stumian,  i.e. 

the  number  of  nodal  points  increases  or  decreases  mono- 

tonically  as  oj  ■*  oj   .   It  will  be  shown  in  chapter  V  that 

o 

this  property  is  lost  in  the  steady  flow  equilibrium. 

The  tip  of  the  cusp  continuum  can  be  treated  in  an 
identical  way  and  yields  similar  results.   We  only  have  to 
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observe  that  det  A  vanishes  when  C  vanishes,  and  that 
although  the  inner  region  equation  is  not  the  Legendre 
equation,  it  still  depends  on  z   alone. 
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IV.  3   A  Necessary  and  Sufficient  Criterion  for  Local  Stability 

We  have  established  that  only  one  point  of  the 
continuous  spectrum  can  serve  as  an  accumulation  point 
of  complex  eigenvalues,  and  this  is  the  overlapping 
point  of  the  Alfven  and  cusp  continua.   To  achieve 
local  stability,  i.e.  non  existence  of  such  a  sequence, 
one  of  the  following  three  conditions  must  be  satisfied: 
(i)    The  (finite)  singularities  of  equation  ( c) /  for 
complex  co  near  to  (the  point  of  overlap)  are  in 
one  side  of  the  real  axis  in  the  x  plane.   i.e. 


(42)  h2  -  pa2  <  0 


(ii)   A  non-oscillation  condition,  requiring  that  when 

we  set  uj  =  co  ,  the  real  solutions  of  (C)  do  not 
o 

oscillate  near  x  =  0,  i.e. 


(43)  a  <  I 


with  a  defined  in  (8) . 

(iii)   If  the  two  cusp  singularities  are  in  the  same  half 

2      2 
plane,  upper  or  lower,  i.e.  Bh   <  pa    ,  then 

either  (31. b)  or  (31. c)  are  satisfied.   If  (31. b) 

holds,  there  will  exist  an  accumulating  sequence 

for  some  location  of  the  boundary  and  not  for  others 
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We  now  write  these  conditions  in  terms  of  the 
plasma  profiles.   We  recall  that  the  coefficients  (10) 
of  the  equation  were  evaluated  at  the  point  were 

—  B„+kB   =  0,  which  determines  the  ratio  r  in  terms  of 

r   8    z  k 

the  magnetic  field.   We  further  notice  that  the  coef- 
ficients (10)  are  homogeneous  of   degree  zero  in  m  and 
k,  and  thus  are  functions  of  j-  only.   It  is  therefore 
possible  to  eliminate  the  wave  numbers  from  the  equation, 
so  that  all  the  stability  conditions   can  be  expressed 
in  terms  of  the  plasma  profiles  only. 

Be 

We  define  y  =  —5—  which  is  essentially  the  pitch 
rts 
z 

of  the  magnetic  field.   y'  is  then  the  shear  of  the 
magnetic  field,  and  after  expressing 


h  =  -kB   ^— 

z  y 


a  =  -k(i  Q'    -   w') 


condition  (42)  takes  the  form 


\2    u'2  <o(Q'    -    uw1)2 


(44)  By'   <p(ft'  -  yw') 

z 


with  0.   =   —  and  w  =  u   .   The  right  hand  side  of  (44)  is 
r  z 

essentially  determined  by  the  shear  of  the  equilibrium 
flow,  and  we  see  that  with  large  sheared  flow,  local 
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modes  are  impossible.   In  particular,  every  local 
instability  can  be  stabilized  by  a  sheared  axial  flow 
which,  as  we  recall,  can  be  introduced  into  every 
equilibrium  state  without  changing  the  other  quantities, 
The  stabilizing  effect  of  sheared  flow  can  be  explained 
by  the  plasma  particles  being  carried  away  from  each 
other,  thus  preventing  localized  modes  from  occurring. 
We  define  the  non  dimensional  quantities 


V 


Be 


Inequality  (43)  takes  the  form 
^rB^  (iVxlL')2  -  W\   I'    +  \r2   ^  (pa*,.  +  p. 

Be 

(45)  R2 

>   — ^  ~   [(1-B)  [j  4>A{l-<p2)-4>3<p+  j4>2ip2]+B[^p2(l+<p2)-2M+$2)] 

When   we    let    the    flow   vanish,     (45)    reduces    to 


(46)  ^   rB^     (-^V    +    p'     >    0 


which  is  Suydam's  condition  [10].   Condition  [44]  can  be 

2 
written  as  $   >  1  and  therefore  (45)  may  be  considered 

2  2         2 

only  for  <p      <  1.      $   =  1  or  $   =  6,  our  analysis  fails 
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and  the  singularity  at  infinity  of  the  inner  region 
equation,  becomes  an  irregular  singularity.   The 
analysis  also  fails  in  the  case  of  a  zero  shear 
magnetic  field  (p*  =  0),  which  requires  a  separate 
treatment. 

An  interesting  case  is  when  the  flow  is  purely 
axial.   i.e.,  ty   =  ft  =  0.   In  this  case  (45)  reduces  to 

(47)       ±rB*  ,  ^,2  ,1-;2,  +  p.  >  2-«4     4 

When  the  shear  of  the  flow  is  small  enough  so  that 
B  -  y   >  0  ,  (47)  is  more  difficult  to  satisfy,  compared 
with  (46) .   Thus  the  introduction  of  flow  into  the 
equilibrium  can  cause  local  instability,  although  large 
flow  shear  makes  the  plasma  locally  stable. 
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CHAPTER  V:   The  Lagrangian  Description 
V.l   The  Equations  of  Motion 

The  Lagrangian  formulation  of  the  ideal  MHD  equations 
"near"  a  steady  state  equilibrium,  was  carried  out  by 
Frieman  and  Rotenberg  [13]  and,  as  will  become  apparent 
shortly,  turns   out  to  be  superior  to  the  Eulerian  descrip- 
tion wherever  stability  analysis  is  involved. 

The  position  vector  r  of  a  plasma  particle  which 

was  at  r   at  t  =  0 ,  is  written  as 
—  o 


(49)  r  =  r°  +  £(r°,t) 


where  r   describes  the  equilibrium  trajectory  of  the 

same  particle  (located  at  r   initially),  and  E,  (r  ,t) 

— o  —  — 

is  the  displacement  from  equilibrium  -  assumed  to  be 
small.   The  independent  variables  are  taken  to  be 
r  ,  t  rather  than  r  , t  so  that  the  equilibrium  quantities 
are  time  independent.   After  linearization,  the  conser- 
vation of  mass,  entropy  and  magnetic  flux,  yield  to 
first  order  in  £: 


p(r°  +  i)  =  p(r°) (1-  V°.£) 


p(r°  +  £)  -  p(r°)  -  pa2  V°-  £_ 


B(r°  +  U  =  B(r°)  -  B  V°.£  +  b-V°  £ 
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where  V   is  the  gradient  operator  in  the  r   coordinates 
The  displacement  vector  £  satisfies  the  equation 

(50)  p£tt+  2A^t  +  HI  =  ° 

where 

AC  =  pu' V£ 


■H4  =   V(pa2V-C  +  C"Vp-B-Q)+B-V  Q 


+Q- V  B  +  V-  (PC  u-V  u  -  pu  u-V  K) 

Q  =  curl(£  x  b) 

All  coefficients  are  the  unperturbed  quantities  and  are 
functions  of  r   only.   The  superscript  was  dropped  every- 
where for  simplicity. 

Equation  (50)  is  an  initial  value  problem  and  the 
proper  values  to  prescribe  are  C.fC..  at  t  =  0 .   A  well 
posed  initial  value  problem  for  the  Ludquist  equations 
(L)  requires,  that  we  give  the  values  of  p,p,u,B 
initially.   It  is  clear  even  from  counting  the  number 
of  data,  that  the  Eulerian  and  Lagrangian  descriptions 
are  not  entirely  equivalent.   It  was  pointed  out  to  the 
author  by  H.  Grad,  that  the  class  of  flows  resulting  from 
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prescribing  the  Lagrangian   data,  is  contained  (in 
the  strong  sense)  in  the  Eulerian  flows.   An  illustra- 
tive example  is  the  "frozen  in"  field  equation 
B   +  curl(B  x  u)  =0,  where  B   is  a  static  equilibrium 
field,  and  B,  u  are  the  perturbations.   A  change  to 
Lagrangian  variables  requires  that  we  write  u  =  £. 
and  integrate  the  equation  with  respect  to  time. 
The  result:  B+  curl(B   x  £)  =  [b  +  curl  (B<^x£)  1  t=0  • 
The  right  hand  side  is  always  dropped  from  the  equation, 
but  this  imposes  a  (non-holonomic)  constraint  on  the 
flow,  namely  B(t=0)  =  curl(£(t=0)  x  B  ).   An  important 
outcome  of  this  constraint  is  that  an  unperturbed 
system  with  closed  line  magnetic  field,  will  remain  so 
under  the  Lagrangian  perturbation,  while  the  Eulerian 
data  allow  the  freedom  of  changing  the  topological 
character  of  the  field  lines. 

We  return  to  the  question  of  stability  of  equation 
(50).   By  choice  of  the  independent  variables,  time  does 
not  appear  explicitly  in  (50)  and  the  equation  can  be 
Laplace  transformed  in  time,  or  equivalently  we  can  as- 
sume an  eluJ  dependence.   Thus,  (50)  takes  the  form 


(51)  -iu2p£  +  2wiA  4  +  H£  =  0 
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A  natural  Hilbert  space  setting  for  this  problem 
can  be  obtained  by  introducing  the  scalar  product 


U,H)  = 


-,*  . 


i   dV 


with  the  integration  carried  over  the  volume  of  the 
plasma,  and  n_*is  the  transposed  and  complex  conjugate 
vector  of  n_  .   It  can  be  easily  shown  that  for  every 
K    whose  normal  component  vanishes  at  the  boundary  of 
the  plasma,  the  quadratic  forms  (H£,£)  and  (iA£,  £) 
are  real.   This  proves  that  the  operators  H  and  iA 
are  symmetric  (formally  self  adjoint) .   We  can  see 
now  how  non  self  adjointness  sets  in  when  flow  is 
introduced  into  the  equilibrium  state.   Without 

flow,  A  =  0  and  (51)  takes  the  form  of  a  classical 

2 
eigenvalue  problem  with  the  result  that  lo  must  be 

real.   In  the  presence  of  flow,  no  such  simple  result 

is  available,  although  regions  in  the  complex  plane 

free  of  "eigenvalues"  u ,  can  be  found  without  much 

labor.   There  is  a  considerable  amount  of  literature 

dealing  with  the  non-linear  eigenvalue  problem  (51) , 

and  an  extensive  reference  list  can  be  found  in  reference 

14 .   Most  of  the  results  deal  with  variational  properties 

of  the  eigenvalues,  completeness  of  the  eigenvectors,  etc., 

but  always  under  restrictive  assumptions  which  assure 

the  reality  of  to  (e.g.  H  is  positive),  or  deal  with  finite 

dimensional  spaces. 
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The  stability  question  received  less  attention, 
but  from  what  is  known   it  seems  highly  unlikely  that 
a  simple  criterion  for  stability  can  be  found  for  the 
general  case.   As  an  illustration  to  the  difficulty 
encountered,  we  state  the  condition  derived  by  Barston 
[15]  for  the  stability  (boundedness  in  time)  of  the 
system  u    +  Au   +  Hu  =  0  in  an  n-dimensional  Eucledian 
space,  with  H  and  iA  Hermitian,  constant  matrices. 
The  system  is  stable  if  and  only  if  there  exists  a 
positive  definite  Hermitian  matrix  p  and  an  anti- 

Hermitian  matrix  a,  such  that  A  =  pa  +  ap  , 

2    1 
H  =  p(a   +  t")P-   Although  such  a  condition  is  useful 

for  constructing  stable  systems,  it  is  practically 

useless  as  a  stability  test. 

We  continue  now  with  the  special  case  of  the 

screw-pinch  configuration.   Initially,  the  equilibrium 

profiles  are  assumed  to  depend  on  r   alone.   Since  the 

radial  velocity  vanishes,  we  have  r   =  r   and  the  co- 

efficients  of  (51)  do  not  depend  on  0  , z  .   We  Fourier 

transform  the  equations  with  respect  to  these  variables, 

and  get 


£.(r°f  t)  =  £(*0)  exp  i  (m0°+kz°+u)t) 


We  can  express  8   ■  0   +       ,z=z+wt,  and  in 
^  o  o 
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terms  of  the  initial  location  of  the  fluid  particle, 
we  have 


£(r°,t)  =  £(rQ)  exp  i(m8Q  +  kzo)e10t 


where  as  before  a  =  uj+k-u  .  According  to  the  usual 
Lagrangian  procedure,  the  velocity  field  is  obtained 
by  differentiating  r  with  respect  to  time,  keeping  r 
fixed.   By  differentiating  (49)  we  get 


u(r,t)  =  uo(r°)+(uo.V°)£  +  £t 


£q  fixed 


here  u   is  the  equilibrium  flow  and  u  is  the 
total  flow  at  r  (not  r  ) .   By  a  Taylor  expansion  we  have 
u(r)  =  u(r°)  +  (£_•  V°)u   to  first  order,  and  therefore 


(52)         6u  =  u(r°)-u  (r°)  =  ia£-£r(rB,') 


where  ftr  ,  w  are  the  angular  and  axial  equilibrium  flow. 
6u  is  the  velocity  perturbation,  previously  called  u  , 
and  (52)  provides  the  transition  from  Lagrangian  to 
Eulerian  analysis.   In  particular,  the  radial  components 
are  related  simply  by  6u  =  ia£  ,  which  was  used  in 
the  derivation  of  equation  (B) .   The  conclusion  is 
that  (51)  can  be  reduced  to  a  single  second  order 
differential  equation  in  £   ,  which  is  exactly  equation  (B) 
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In  view  of  future  benefits  we  will  nevertheless  re- 
derive  this  equation,  using  the  Lagrangian  formulation 
only . 

We  rewrite  equation  (51)  in  the  form 


(53)  H  £  =  0 

co — 


2 

where  H   =  -pco   +  2coiA  +  H 

co 


claim:  For  co  real,  equation  (53)  has  the  equivalent 
variational  formulation 


6(H  £,£)  =  0 

co 


where  the  variation  is  with  respect  to  all  twice 
dif f erentiable  complex  vectors  £_,  with  a  vanishing 
normal  component  on  the  boundary  of  the  plasma, 
proof:   For  co  real,  H   is  a  symmetric  operator. 


6(H  CO  =  (H  6£,0  +  (H  £,60  =  2  Re  (H  £,«€>  =  0 

10  CO  CO  CO 


By  taking   i6£  instead  of  6E,   we  find  Im(H  £,6£)  =  0. 

co 

Thus  H  E,    =  0. 
co- 
in the  screw  pinch  case  (H  £,£)  contains  r-deriva- 

h) 

fives  of  ■'.   only,  and  is  algebraic  in  f,Q,E,z.      Moreover, 
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if  we  let 


C  =  Cr  ,  n  =  Ue   ,5  =  Hz 


then  the  coefficients  of  the  quadratic  form  are  real, 
and  (53)  has  a  real  solution.  We  thus  consider  £,n,C 
to  be  real  and  view  (H  £,£)  as  an  algebraic  quadratic 

CO 

in  n t C  •   This  quadratic  has  the  following  terms: 
Second  degree  terms  in  n ,  t, : 


(54)   s(n,c)  =  pa2(^n  +kc)2+K2(nBz-;BQ)2  -po2(n2+c2) 


First  degree  terms: 


|(jh+kC)  [pa2(r£)'    +    pv2a+2(nBz-cBe)[p(kBe+   ^Bz)     +  (^-kBg)  £  '  ] 


-4pft0£ri 


Zeroeth  degree  terms : 


B2  C'2  +  ?2[H2-pa2  -(£f_)'  + 


r,   2        „   2 

B       B„ 


r      r 


This  quadratic  is  integrated  along  the  range  of  the 
plasma  with  the  measure  rdr .   The  terms  containing  ££' 
were  eliminated  by  integration  by  parts.   The  terms 
containing  r\ ,  c,   can  be  written  as  I  (y)  =  (Sy  ,y) +2  (b,y) 
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where  y  =  (";)  ,  S  is  the  matrix  associated  with  (54), 
and  b  is  a  vector  linear  in  £,£'.   I (y)  is  an  algebraic 
quadratic  and  it  is  stationary  for  y  ,  where 


(55)  y   =  -S  1b   . 

o 


Substituting  these  values  of  n / C  we  get 


(V'e)  = 


'2      2 

(PC  *    -   qC  )  dr 


and  the  corresponding  Euler-Lagrange  equation  is 

(pO  *  +  qc  =  o  . 

The  coefficients  p,q  are  exactly  those  defined  in  (B) . 
The  substitution  (55)  was  instrumental  in  reducing  the 
original  equation  (53)  to  the  scalar  equation  (B)  in  £ 
alone.   Although  the  reduction  was  based  on  the  assump- 
tion that  w  was  real,  we  know  that  the  coefficients  of 
(53)  are  analytic  in  w,  and  thus  we  can  conclude  that 
equation  (B)  applies  to  cu  complex  too,  as  is  known  from 
the  Eulerian  description. 
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The  matrix  S  has  the  determinant 


2  4    2 
det  S  =  p  c  +CK   =  A 


as  defined  in  (B),  and  it  is  now  clear  how  this  apparent 
singularity  enters  into  (B) .   From  the  discussion  in 
section  II.  2  we  know  already  that  E,   has  no  singularity 
at  the  point  where  A   =  0.   As  to  n. ,  z,  >    we  use  (52)  to 
express  them  in  terms  of  the  Eulerian  variables  which 
are  not  singular  at  this  point,  and  we  conclude  that 
A, 2  is  not  part  of  the  spectrum  of  the  Lagrangian 
equation,  despite  (55) .   We  also  find  that  a   =  0  is  not 
part  of  the  spectrum  in  this  formulation.   This  dif- 
ference between  the  Eulerian  and  Lagrangian  spectra  is 
not  surprising  in  view  of  the  discussion  at  the  begin- 
ning of  the  chapter.   The  continuous  spectrum  is  related 
to  solving  an  initial  value  problem,  and  this  problem 
is  not  equivalent  in  the  two  descriptions. 
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V.2   Properties  of  the  Point  Spectrum 

We  take  the  inner  product  of  (51)  with  £  and  solve 
for  u>. 


1 


(56)  w  =  T^TT   <<«'«*  /(±a5,€)2+<h5,5HpC.5>] 


One  observes  [13]  that  a  sufficient  condition  for  stability 
(a)  real)  ,  is 


(57)         (iACU2  +  (H€,£)  (pC,€)  >  0 


for  all  E, .      This   is  known  as  the  overdamping  condition , 
If  oj  is  a  complex  eigenvalue,  we  have 


(58)  Re  u  =  UAg,£) 

(P£,C) 


(59)  m2  =  zm^L 

(pC/C) 


where  £  is  the  corresponding  eigenf unction.   These  results 
were  derived  by  Bars ton  [16] .   It  was  also  observed  by  the 
same   author  that  if  H  is  bounded  below,  say  (HE,,E,)>_   Mp£»5) 
for  all  ^(understood  to  have  a  vanishing  normal  component  on 
the  boundary)  and  with  X    >    -<*>    ,  we  have 
(60)  lu)l2  <  -  X 
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This  condition  restricts  the  complex  eigenvalue  to  a 
circle  in  the  complex  plane,  and  what  is  more  important 
from  the  point  of  view  of  stability,  it  provides  an 
upper  bound  for  the  growth  rate  of  instabilities.   A 
better  bound  for  the  growth  rates  can  be  derived  when  A 
is  a  bounded  operator  (as  will  be  shown  to  be  the  case 
in  the  screw  pinch  equilibrium) .   We  make  a  real  trans- 
lation of  the  complex  u>  plane 

u>  =   a   +   u ' 


with  a  real  constant  a.   We  get 


H   =  -pa)'   +  2co'  (iA-ap)+H 
a)  a 


where  Ha  is  defined  in  (53) .   The  last  expression  is 
similar  in  form  to  Hi,  and  using  an  argument  similar 
to  (60)  we  can  conclude 


(61)  lor  I2   <  -   A 

i   i    _      a 


a   =   Inf   


=   Inf 

S    <P5,€) 


One  can  pick  an  optimal  a  by  maximizing  X    .   We  then  get 

2  (HttS»5) 

(62)      (Im  cj)z  <  -  Max     Inf    — 


a  real    £     (p£#S) 
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for  any  complex  eigenvalue  of  (53) .   Our  new  bound 
enables  us  to  state  the  following 

Corollary:   A  sufficient  condition  for  stability  is 
that  the  following  holds: 


(H  C,C) 
(63)  Max    Inf   ,  -  rN   >  0 

a  real   ?    {^'° 


The  next  section  will  be  devoted  to  the  interpretation 
and  investigation  of  this  criterion.   We  now  show  that 
in  the  screw  pinch  equilibrium,  the  last  criterion  is 
applicable (i .e.  the  previous  assumptions  are  satisfied) 
Claim :   A  is  bounded,  H  is  bounded  below. 
Proof:   BY  the  definition  of  A 


(64)       AC  =  pu-V  £  =  ip(k-u)£  +  pQ(    £r) 


A  is  not  a  differential  operator,  and   |a|   <_  Max  p[|k-u|  +  |ft|] 
Changing  our  notation  from  a  to  aj(real)  ,  we  now  want  to 
show  that  H   is  bounded  below.   The  claim  refers  to  the 
special  case  w  =  0.   From  the  explicit  expression  of 

,i)    detailed  in  section  V.l,  we  see  that  this  quadratic 
can  be  viewed  as  an  algebraic  quadratic  in  " ' ,  having  as 
first  and  second  degree  terms  the  expression 


(65)     JU')  =  B*  C,2  +  2f-,'[pa2(pn+kO  +  (kBn-  p-B^  (B()  .'.-Bzn)  1 
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All  other  terms  are  of  second  degree  in  £  ,  n  »  C  •   The 
coefficients  of  the  quadratic  are  real  and  minimization 
can  be  carried  out  by  considering  only  real  vectors 
(£,n#0.   The  reason:  if  (H  £,£)   >  A  (p£,0  for  all 

CO  —        GO 


real    £ ,    then 


(H    (C+i<|)),C+i<t»)    =    (H   £,£)+(H   d>,<f>)  -2    im(H   E,<|>)    =A    [  (p?,£)  +  (p<p  ,$) 

CO  CO  00  CO  GO 


because  Im(H  £,,$)    =  0  when  £  and  £  are  real.   Similarly 
(p(£+i<f>),  £+i<}>)  =  (pC/S)  +  (p<{>#<J>)  t    and  we  see  that  the 
same  bound  applies  to  real  and  complex  vectors.   We 
therefore  assume  C/H/C  to  be  real,  and  minimize  (65) 
with  respect  to  E, '  .      The  minimum  value  of  J(^')  is 

(66)    -  -|  [pa2(^  n+kc;)  +  (kBe-  ^  (BqC-B^)]2  . 
B* 

2     9 
(It  is  assumed  here,  as  in  section  II. 2,  that  B^ ,  pa  do 

not  vanish  in  the  range  of  the  plasma) .   To  find  a  lower 

bound  for  (H  £,£)  we  replace  J(C')  by  (66) ,  and  get  an 

GO 

algebraic  quadratic  in  E,,r\,c,.      Such  a  quadratic  is  bounded 
below  by  the  lowest  eigenvalue  of  the  associated  matrix 
and  the  claim  is  proved. 

We  carry  out  this  minimization  in  order  to  obtain  a 
lower  bound  for  A  ,  and  thus  an  upper  bound  for  the  growth 
rates.   After  simple  algebraic  manipulations  we  write  the 
expression  as  a  sum  of  two  quadratic  forms,  the  first  of 
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which  can  be  represented  by  the  symmetric  matrix  Q 
with  elements 


2 


Be      B*  /   h2 

Q    =  [r(-^y)  '  -  (— -)  '  +  H^]/[l  +  /l-  -4-=-  ] 
r         r  kb 


/      ?    l 
Q12  =  p(pa2  +  PV2)/[l  +  A  -  -4J-J  ]2 

K  B 

(67)       R      m         /    ^2"  I 
Qi3  "  F(kBe  +  FBZ)/[1  +  A  "  J^]2 

K  B 


Q22  =  6B2K2 


Q23  =  6BK(kBe  -  ^Bz) 


Q33  =  (1-6)H2+6B2K2 


,  x 

where  Q  is  represented  in  the  space  of  vectors  x  =  (y) 

—  z 

with 


x  =  [1  +  /1-  -4 


2    2, 

1  /m   ...        1 


(68)   x  =  [1+/1-  -J-j  ]  £   ,  y  =  i(|n  +kO  »  z  =  t(nBz-CBQ) 

K  B 


We  have 


H   ,U  >  (Qx,x)  -p[o2(^2  +  n2  +  c2)+4no?n] 


The  minimum  of  the  right  hand  side  would  not  increase  if 
we  minimize  the  two  terms  separately.   From  (68)  we  find 


(1-/-  «    (n2n2)<  y2  +  22  <  (1+/-  -|-5  (n2+C2) 
ITB*  KB 
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(which  is  responsible  for  the  definition  of  x) .   Thus, 
if  q  is  the  lowest  eigenvalue  of  Q  (and  we  assume  it  to 
be  negative  at  least  for  same  r  in  the  plasma) ,  we  find 

(69)     Inf  ,  "     I   Min  [SL(i+/i-_"   )  -  (o%2  |  a  |  |  ft  |  )  ] 
£  KQK'V  r    P        KZBZ 

If  q  is  positive  everywhere,  we  replace  the  radical 
term  by  its  negative  value.   The  minimum  over  r  is  taken, 
because  a  vector  £_  can  be  chosen  to  be  concentrated 
around  the  point  where  the  minimum  is  attained.   Expression 
(69)  has  to  be  optimized  with  respect  to  all  real  u>. 

It  is  possible  to  save  much  effort  by  choosing  cj  =  co 

2 
such  that  Max  (a  +2  |  a  |  |  ft  |  )  be  minimized  by  to  . 

Relation  (69)  provides  us  with  the  following  result. 

If  K  is  large  enough  and  H  =(=  0  in  the  range  of  the  plasma, 

then  the  matrix  Q  is  positive  definite  and  we  get 


(70)        (Im  a))2  <  Min  Max(a2+2|a|  |ft|  ) 


for  every  complex  eigenvalue  00. 

Before  continuing  the  discussion  about  the  best  bound 
X      1    we  find  the  range  of  Re  00  for  complex  eigenvalues. 
From  (64)  we  have 


(iA£,5)  =  - 


p[(k-u) (£2+n2+C2)+2ft£n]rdr 
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The  eigenvalues  of  the  quadratic  under  the  integral  sign 
are  -k'U  ,  -k*u±fl  ,  and  by  (58)  we  get 


(71) 


Min[-k-u    -    \Sl\]    <_  Re   oi   <   Max    [-k-u   +    |Q|] 


w 


As  a  demonstration  we  use  the  results  just  established, 
for  the  special  case  of  solid  rotation,  i.e.  Q   =   constant, 

=  0.   In  this  case  o  =  w+mft  and  the  optimal  oj  for 
inequality  (69)  is  u>   =  -mft,  so  that  o  (w    )  =0.   Relation 
(70)  yields  stability  for  large  K  (and  H=j=0)  .   In  the 
case  of  instability  we  find  from  (71) 

-  |  ft  |  £  Re  w  +  mft  <_  |  ft  | 

We  return  now  to  the  optimization  (with  respect  to 
w)  of  A  .   If  H   was  a  self  adjoint  operator  (and  not  only 
formally  self  adjoint)  ,  we  could  have  concluded  that  A 
was  the  lowest  point  in  the  (A)  spectrum  of  H  ,  character- 
ized by  the  equation 

(72)  H^   =  Xp£ 


Nevertheless,  our  purpose  is  to  suggest  a  candidate  oj 
for  optimization  and  not  to  actually  find  X  ,  and  to 
this  end  we  can  treat  H   as  a  self  adjoint  operator. 
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The  continuous  spectrum  of  H   is  easily  detected. 
We  repeat  the  procedure  carried  out  at  the  beginning 

of  the  chapter  in  reducing  (53)  to  (B) .   The  only  dif- 

2  2 

ference  is  that  go   is  replaced  by  w  +A .   We  thus  reduce 

(72)  to  a  second  order  ODE  similar  to  (B) ,  and  the 

continuous  spectrum  exposes  itself  through  the  vanishing 

of 


-22        ~      22 

A  =  H  -pa  -PA  ,   c  =  3H  -pa  -PA 


We  find  that  the  continuous  spectrum  is  composed  of  all 
points 


a    1  u2     2     .    6  „2     2 
A  =  —  H   -a   or  A  =  —  H   -a 

P  P 


for  some  r  in  the  plasma  .  The  lowest  point  of  the 

6   2   2 
continuous  spectrum  is  clearly  A  (oo)  =  min  (—  H  -a  )  . 

c       r   P 

If  there  is  not  any  eigenvalue  lower  than  A   (for  all  go)  , 

we  can  choose  go   to  be  the  go  maximizing  A  .   We  now  assume 

o  ^   c 

that  there  exists  an  eigenvalue  A   <  A   .   From  (B)  we 
know  that  every  eigenvalue  A  is  an  analytic  function  of 
go  as  long  as  A  (go)  does  not  cross  the  continuous  A-spectrum. 
Another  observation  we  make  is,  that  the  eigenvalue 
curves  A  (go)  do  not  cross  each  other.   The  reason:   at  a 
crossing  point  there  will  be  two  eigenf unctions  for  one 
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A,  which  is  impossible  because  A  is  determined  by  a 

second  order  differential  equation.   As  a  result  we 

can  conclude  that  the  lowest  eigenvalue  A  (oj)  is  an 

o 

analytic  function  of  oj.   We  can  now  differentiate 
(72)  with  respect  to  oj. 


(73)  H  £  +  H  £  -  X  p£  +  A  PC 

00 —     CO —     o  —     o  — 


where  x  =  y—   .   Multiplying  (73)  by  £  ,  and  using  the 

a  OJ  — 

fact  that  £  satisfies  the  homogeneous  boundary  conditions 
(B1-B2) ,  we  get 


(74)  (2u)+AQ)  (p£,£)  =  2(iA£f£) 


We  are  interested  in  the  optimal  A   =  A  (co  )  ,  and 

o    o   o 

we  find  from  (7  4) 


(iA£,£) 
wo  '   (p£, O 


Similarly  to  (71) ,  we  get 


(75)       Min[-k-u  -  Iftl]  <  id   <  Max  [-k-u  +  Ini] 
r    —  —        —   o  —  r    —  — 


which  is  expected  in  view  of  (71) .   We  want  to  show  that 

if  \   =  A   ,  then  there  exists  a  unique  optimal  a>  . 
0     0)  ^  o 
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It  would  be  sufficient  to  show  that   A   <  0.   Dif- 

o 

ferentiating  (73)  aga:.n,  we  get 


HE  +  2H  E+C  =  A  o£  +  2pA  E  +  A  pE 

co —      co —  —    o  —      o —    o  — 


Multiplying  by  E_  ,  and  using  the  relation  H   =  -2p  , 
we  find 


4((iA-cop)E,E)  -  2(pE,E)  =  AQ(pS,E)+2Ao  (pE,E) 


We  now  multiply  (73)  by  2E_  and  subtract  it  from  the 
last  equation 


A  (p£,E)  =  -  2(pE,E)  +  2A  (pE,E)-2(H  E,E) 

O  O  CO 


If  indeed  A   =  A   ,  then  the  last  two  terms  give  a 

o      CO 

negative  contribution,  and  we  conclude  that  A   <  -2 . 
The  optimizing  co   is  thus  unique. 

The  next  section  will  offer  a  simple  way  to  compute 

the  lowest  eigenvalue  of  H   numerically.   The  analytical 

co 

results   we  derived  here  can  thus  have  a  practical  use. 
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V.3   A  Sufficient  Condition  for  Stability 

In  the  previous  section  we  derived  a  sufficient 
condition  for  stability: 


Treal   T   (pC  '° 


(63)  Max      Inf    "r  _  .   >  0 


2 
In  the  static  equilibrium  case  H   =  H-pco  ,  and  the 

maximum  with  respect  to  cu  is  obtained  for  oj  =  0 .   The 
resulting  sufficient  condition  for  stability  is  H  >  0  , 
which  is  known  to  be  a  necessary  condition  as  well. 
This  is  not  the  case  when  flow  is  present. 

Condition  (63)  was  originally  derived  (using  a 
different  method  than  ours)  by  Bussard  [17] ,  and  was 
believed  by  him  to  be  a  necessary  and  sufficient  con- 
dition for  stability.   The  proof  given  by  this  author 
applied  only  to  finite  dimensional  spaces,  and  of  course 
the  proof  of  the  necessity  part  was  erroneous.   Bussard 's 
idea  however,  offers  a  visual  interpretation  of  condition 
(63)  as  well  as  a  demonstration  of  its  unnecessity.   Let 
us  then  consider  the  case  p  =  1  and  iA,H  symmetric  constant 

matrices  of  order  n.   The  Eigenvalue  equation  H  t,   =  0 

Co 

has  now  exactly  2n  solutions,  given  by  det(H  )  =  0.   The 
system  would  be  stable  if  we  can  demonstrate  the  existence 
of  2n  real  eigenvalues.   To  this  end  we  consider  the 
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A-eigenvalue  problem 


H^  =  A  (U)e 


By  standard  theorems  [18]  we  know,  that  the  n  eigenvalues 
can  be  described  by  n  functions  A,  (to)  (k=l,...n)  which 
are  analytic  in  a  neighborhood  of  the  real  line.   Starting 
at  co  =  0  ,  for  example,  we  label  the  eigenvalues  *k(0) 
according  to  the  natural  order, with  A.,  the  smallest  one. 
As  we  let  a)  vary  along  the  real  axis,  a  curve  ^v(^)  may 
cross  the  horizontal  axis  (in  the  (go, A)  plane).   The 
point  of  intersection  is  an  co-eigenvalue  because  the 
equation  H  £  =  0  has  a  non -trivial  solution.   All  real 
eigenvalues  co  are  given  by  such  points  of  intersection. 
The  curves  A,  (co)  behave  like  -co   for  |  co  |  large  and  thus 
are  negative.   If  for  some  co  A,  becomes  positive,  there 
will  be  at  least  two  points  of  intersection.   If  all  of 
the  curves  become  positive (fig .  8) 
we  clearly  have  2n  real 


co-eigenvalues  and  the 

system  is  stable. 

Condition  (63) 

states  that  if  the  pI( 

lowest  A-eigenvalue  becomes  positive  for  same  co,  we  have 

stability.   Two  counter  examples  would  show  why  this 


>  co 
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>U) 


>  u 


FIG.  9  FIG.  10 

condition  is  not  necessary.   In  (fig.  9)  the  eigenvalue 
curves  intersect.   The  lowest  eigenvalue  for  any  given 
id  is  always  negative,  while  the  analytic  curves  A^do) 
intersect  the  w  axis.   Parenthetically  we  add  that 
in  the  screw  pinch  case  this  situation  does  not  occur. 
The  second  counter  example  is  when  one  of  the  higher 
curves  crosses  the  axis  more  than  twice  (fig.  10), 
where  we  can  still  have  stability  although  A,  (uj)  is 
always  negative. 


We  return  now  to  the  minimization  of 


<hbcc> 


(pCO 

for  fixed  u.   From  the  previous  section  we  know  that 

there  is  hope  to  find  H   non-negative  only  if 

2      2 
BH   -  po   >  0  for  all  r  in  the  plasma.   We  derive  the 


same  result  in  a  more  direct  way  by  actually  minimizing 
(H  CO 

—7 — b — p-r-      .   In  fact,  because  we  are  only  interested  in 
U»CO 

the  non-negativity  of  H   ,  we  can  take  the  internum  of 


ui 


(H  £,£)  without  normalizing 


U 


Like  (54),  (H  f. ,-".)  is  an  algebraic  quadratic  in 


u 


n,C  and  we  first  minimize  it  with  respect  to  these  two 
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variables.   The  minimum  is  finite  if  and  only  if  the 
second  degree  term  S(n»0  is  positive  definite,  or 
if  it  is  semi  definite,  then  the  first  degree  term 
also  vanishes  when  the  second  degree  term  does  so. 
The  positivity  of  S  requires  that 


(76)  trace  S  =  K2B2  -  2pa2  >  0 


det  S  =  p  a  +K  B*(BH  -po  )>  0 


be  satisfied  for  all  r.   From  the  two  inequalities  it 
follows  that  H2-po2  >  0. 

Assuming  S  to  be  positive  definite,  we  minimize  with 
respect  to  r\ ,  t,    and  get  the  known  form 


f 


(77)  FU)  = 


'2      2 

(PC    -  qC  )  dr 


with  p,q  defined  in  (B) .   F(£)  is  bounded  below  only  if 
p  >  0  everywhere  in  the  plasma.   Rewriting  p  as 


rB,2(H2-pa2) (0H2-pa2) 


P  = 


det  S 


it  is  clear  that  in  addition  to  (76)  we  need  the  condition 


(78)  BH2  -  pa2  >   0 
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to  hold  everywhere.   Inequality  (78)  is  more  restrictive 
than  (76)  (unless  H  =  a  =  0)  and  it  is  sufficient  to 
oo/.sider  only  it.    w  can  still  be  chosen  arbitrarily  as 
long  as  (78)  holds,  but  if  H  =  0  at  some  point  r   ,  we 

must  choose  w  =  u  ,  so  that  a(r  ,  co  )  =  0.   The  point 

o  o   o  r 

(r  ,  a)  )  is  recognized  as  the  point  of  overlap  of  the 
o   o 

Alfven  and  cusp  continua.   If  there  are  two  points  in 
the  plasma  where  H  =  0 ,  we  cannot  use  the  sufficient 
condition  unless  a  vanishes  at  both  points  for  the  same 
u>.   This  will  be  the  case  for  instance,  if  the  equili- 
brium velocity  field  is  parallel  to  the  magnetic  field, 
up   to  a  solid  rotation. 

Before  we  continue,  a  certain  difficulty  has  to  be 
removed.   If  the  point  of  overlap  of  the  continua  has 
to  be  chosen  in  order  to  implement  the  sufficient  con- 
dition, the  strict  condition  det  S  >  0  is  violated  at 
that  point.   Luckily  however,  the  first  degree  term  in 
(54)  also  vanishes  at  this  point  and  minimization  with 
respect  to  n ,  C  does  not  yield  -°°  . 

The  form  F(£) (p^O)  in  (77)  is  non-negative  for  all 
E,    satisfying  the  homogeneous  boundary  conditions  (B1-B2), 
if  and  only  if  the  following  condition  holds: 
i)   If  p  is  strictly  positive  everywhere,  then  no  solution 
of  the  Euler-Lagrange  equation  (B)  associated  with  F(£) 
vanishes  twice  in  the  domain  of  the  plasma  [19] . 
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ii)    If  p  vanishes  at  r   corresponding  to  the  overlap 

of  the  continua,  but  otherwise  is  positive  (except 
at  r  =  0) ,  then  no  solution  of  (B)  vanishes  twice 
in  any  of  the  subintervals  (0,r  ),(r  ,R)  [11]. 

It  is  a  simple  task  to  find  numerically  whether  or 
not  this  condition  is  satisfied.   We  can  use  this  cri- 
terion to  compute  A   of  (61)  numerically.   All  we  have 
to  do  is  to  find  the  largest  number  A  such  H-Ap  is  a 

positive  operator.   This  requires  that  we  solve  equation 

2     2 
(B)  with  the  modification  w   ■+  oi  +  A  ,  and  check  if  a 

solution  that  satisfies  the  boundary  condition  on  the 

left  (B2) ,  does  not  vanish  inside  the  plasma.   If  it 

does,  we  take  X  smaller.   The  first  zero  of  the  solution 

in  tne  plasma,   moves  monotonically  (with  A)  toward 

the  boundary  r  =  R,  and  when  it  reaches  it,  A  attains 

the  value  of  X    .   This  procedure  holds  as  long  as  A 

S   2   2 
is  smaller  than  A   =  min(—  h  -a  )  as  described  in  section 

c     r   p 

V.2.   In  this  case  A   is  indeed  the  lowest  eigenvalue 

w    2 
3 
of  H   and  we  know  that  — =-  A   <  0.   This  property  will 

3d) 

be  of  great  help  in  locating  the  unique  optimal  u>. 
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Appendix  B 

All  Eulerian  perturbations  can  be  described  as  linear 
combinations  of  the  three  perturbed  quantities  p^UjU' 
defined  in  section  1.2.   We  write  this  relation  as 

B   =  -u   . 
r   a 

2   2    2 

ipo  (H  -pa  )  -n  =  B<£ 

where 

n.  =  (WBe'Bz'p) 

£  =  (ip* ,u,u' ) 

are  the  perturbed  quantities,  and  B  is  a  5><3  matrix  whose 
entries  depend  on  the  equilibrium  profiles. 

Bll  =  a(Fp°2  "  k2»V 


B12  =  -^B6H"P°2)  l  (H2-P°2)V  -  (H2  +  P°2)^-]  +  J-  kBQH(H2-po2)w' 


+  -  HBA(-BflH-kB  H-PO2) 
r    0  xr  «     z 


B13  =  -B6H(H2-P02) 


2   2 
B_.  =  o(kpo  -K  HB  ) 
1 1  z 
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B22    =    ^HBz[(H2-pa2)v'-(H2+pa2)^]     +    ^-(kBzH-pa2)  (H2-po2)w' 


+  FHBz(FBeH-kBzH  +p°2) 


B23    =    -BzH(H2-pa2) 


B31  =   P°2(F  h_k2V 


B32    =    PBe(H2-po2)  [^-(v1-    p)+kw']     +    |  kBzp2a2v 


-pa(H   -pa    )B'    -   -paHkBQB 
or  u    z 


B33    =    -paBQ(H2-po2) 


B41    =    Pa2(kH-K2B    ) 


mT 


B42         pBzPt (H2~Pa2)v'-(H2+pa2)p]     +kBp(H2-po2)W 


0  0  1  "") 

-pa(H   -po    )B'    +   -paB    (-BQH-kB   H   +    pa    ) 


z        r        z   r 


B 
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•paB    (H2-pa2) 


2       2         2 
B51    =    plT(H    -pa    ) 


B 


52 


p(H2-pa2) [p(^  v'    +  ^  +kw') 
r  r 


-   i-    (-   B0H-kB    H    +    pa2) 
ra      r      9  z 


a '      2 

^  *     -   p'a)] 

a 


B 


53 


o       2         2N2 
—  IH    -pa    ) 
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